diff --git a/.github/dummy_for_action b/.github/dummy_for_action new file mode 100644 index 0000000000..91bdf4d63e --- /dev/null +++ b/.github/dummy_for_action @@ -0,0 +1,2 @@ +#update me to add action comment + diff --git a/.github/jobs/build_docker_image.sh b/.github/jobs/build_docker_image.sh index 5a9b2e8cf7..cefe962fc2 100755 --- a/.github/jobs/build_docker_image.sh +++ b/.github/jobs/build_docker_image.sh @@ -12,3 +12,7 @@ time_command docker build -t ${DOCKERHUB_TAG} \ --build-arg SOURCE_BRANCH \ --build-arg MET_BASE_IMAGE \ -f $DOCKERFILE_PATH ${GITHUB_WORKSPACE} +if [ $? != 0 ]; then + cat ${GITHUB_WORKSPACE}/docker_build.log + exit 1 +fi diff --git a/met/data/config/ConfigConstants b/met/data/config/ConfigConstants index ee0d4855c0..29b2590ada 100644 --- a/met/data/config/ConfigConstants +++ b/met/data/config/ConfigConstants @@ -88,6 +88,7 @@ AW_MEAN = 20; GAUSSIAN = 21; MAXGAUSS = 22; GEOG_MATCH = 23; +HIRA = 24; // Interpolation types NONE = 1; diff --git a/met/data/config/PB2NCConfig_default b/met/data/config/PB2NCConfig_default index 82750e9912..13110d221d 100644 --- a/met/data/config/PB2NCConfig_default +++ b/met/data/config/PB2NCConfig_default @@ -107,7 +107,7 @@ obs_bufr_map = []; // abbreviations in the output. This default map is appended to obs_bufr_map. // This should not typically be overridden. // -obs_prefbufr_map = [ +obs_prepbufr_map = [ { key = "POB"; val = "PRES"; }, { key = "QOB"; val = "SPFH"; }, { key = "TOB"; val = "TMP"; }, @@ -121,7 +121,8 @@ obs_prefbufr_map = [ { key = "D_MIXR"; val = "MIXR"; }, { key = "D_PRMSL"; val = "PRMSL"; }, { key = "D_PBL"; val = "PBL"; }, - { key = "D_CAPE"; val = "CAPE"; } + { key = "D_CAPE"; val = "CAPE"; }, + { key = "D_MLCAPE"; val = "MLCAPE"; } ]; //////////////////////////////////////////////////////////////////////////////// diff --git a/met/docs/Users_Guide/config_options.rst b/met/docs/Users_Guide/config_options.rst index 7f96a62d46..ae5ca26a40 100644 --- a/met/docs/Users_Guide/config_options.rst +++ b/met/docs/Users_Guide/config_options.rst @@ -610,7 +610,7 @@ using the following entries: * MAXGAUSS to compute the maximum value in the neighborhood and apply a Gaussian smoother to the result - The BEST and GEOG_MATCH interpolation options are not valid for regridding. + The BEST, GEOG_MATCH, and HIRA options are not valid for regridding. * The "width" entry specifies a regridding width, when applicable. - width = 4; To regrid using a 4x4 box or circle with diameter 4. @@ -1689,7 +1689,10 @@ This dictionary may include the following entries: * MAXGAUSS for the maximum value followed by a Gaussian smoother * GEOG_MATCH for the nearest grid point where the land/sea mask - and geography criteria are satisfied. + and geography criteria are satisfied + + * HIRA for all neighborhood points to define a spatial + ensemble (only in Ensemble-Stat) The BUDGET, FORCE, GAUSSIAN, and MAXGAUSS methods are not valid for interpolating to point locations. For grid-to-grid comparisons, the @@ -3467,7 +3470,7 @@ of the forecast the observation is used to verify. obs_bufr_map = []; -obs_prefbufr_map +obs_prepbufr_map """""""""""""""" Default mapping for PREPBUFR. Replace input BUFR variable names with GRIB @@ -3478,7 +3481,7 @@ abbreviations to the output. .. code-block:: none - obs_prefbufr_map = [ + obs_prepbufr_map = [ { key = "POB"; val = "PRES"; }, { key = "QOB"; val = "SPFH"; }, { key = "TOB"; val = "TMP"; }, diff --git a/met/docs/Users_Guide/ensemble-stat.rst b/met/docs/Users_Guide/ensemble-stat.rst index ca449cdd9c..f063e23f52 100644 --- a/met/docs/Users_Guide/ensemble-stat.rst +++ b/met/docs/Users_Guide/ensemble-stat.rst @@ -17,6 +17,8 @@ Ensemble forecasts derived from a set of deterministic ensemble members Ensemble forecasts are often created as a set of deterministic forecasts. The ensemble members are rarely used separately. Instead, they can be combined in various ways to produce a forecast. MET can combine the ensemble members into some type of summary forecast according to user specifications. Ensemble means are the most common, and can be paired with the ensemble variance or spread. Maximum, minimum and other summary values are also available, with details in the practical information section. +Typically an ensemble is constructed by selecting a single forecast value from each member for each observation. When the High Resolution Assessment (HiRA) interpolation method is chosen, all of the nearby neighborhood points surrounding each observation from each member are used. Therefore, processing an N-member ensemble using a HiRA neighborhood of size M produces ensemble output with size N*M. This approach fully leverages information from all nearby grid points to evaluate the ensemble quality. + The ensemble relative frequency is the simplest method for turning a set of deterministic forecasts into something resembling a probability forecast. MET will create the ensemble relative frequency as the proportion of ensemble members forecasting some event. For example, if 5 out of 10 ensemble members predict measurable precipitation at a grid location, then the ensemble relative frequency of precipitation will be :math:`5/10=0.5`. If the ensemble relative frequency is calibrated (unlikely) then this could be thought of as a probability of precipitation. The neighborhood ensemble probability (NEP) and neighborhood maximum ensemble probability (NMEP) methods are described in :ref:`Schwartz and Sobash (2017) `. They are an extension of the ensemble relative frequencies described above. The NEP value is computed by averaging the relative frequency of the event within the neighborhood over all ensemble members. The NMEP value is computed as the fraction of ensemble members for which the event is occurring somewhere within the surrounding neighborhood. The NMEP output is typically smoothed using a Gaussian kernel filter. The neighborhood sizes and smoothing options can be customized in the configuration file. @@ -161,6 +163,8 @@ ____________________ The configuration options listed above are common to many MET tools and are described in :numref:`config_options`. +Note that the **HIRA** interpolation method is only supported in Ensemble-Stat. + _____________________ .. code-block:: none diff --git a/met/docs/Users_Guide/reformat_point.rst b/met/docs/Users_Guide/reformat_point.rst index 7a3fce768f..0953afdd7c 100644 --- a/met/docs/Users_Guide/reformat_point.rst +++ b/met/docs/Users_Guide/reformat_point.rst @@ -227,7 +227,7 @@ _____________________ obs_bufr_var = [ 'QOB', 'TOB', 'ZOB', 'UOB', 'VOB' ]; -Each PrepBUFR message will likely contain multiple observation variables. The **obs_bufr_var** variable is used to specify which observation variables should be retained or derived. The variable name comes from BUFR file which includes BUFR table. The following BUFR names may be retained: QOB, TOB, ZOB, UOB, and VOB for specific humidity, temperature, height, and the u and v components of winds. The following BUFR names may be derived: D_DPT, D_WIND, D_RH, D_MIXR, D_PRMSL, D_PBL, and D_CAPE for dew point, wind speed, relative humidity, mixing ratio, pressure reduced to MSL, planetary boundary layer height, and convective available potential energy. This configuration replaces **obs_grib_code**. If the list is empty, all BUFR variables are retained. +Each PrepBUFR message will likely contain multiple observation variables. The **obs_bufr_var** variable is used to specify which observation variables should be retained or derived. The variable name comes from BUFR file which includes BUFR table. The following BUFR names may be retained: QOB, TOB, ZOB, UOB, and VOB for specific humidity, temperature, height, and the u and v components of winds. The following BUFR names may be derived: D_DPT, D_WIND, D_RH, D_MIXR, D_PRMSL, D_PBL, D_CAPE, and D_MLCAPE for dew point, wind speed, relative humidity, mixing ratio, pressure reduced to MSL, planetary boundary layer height, convective available potential energy, and mixed layer convective available potential energy. This configuration replaces **obs_grib_code**. If the list is empty, all BUFR variables are retained. _____________________ @@ -248,6 +248,7 @@ _____________________ { key = 'D_PRMSL'; val = 'PRMSL'; }, { key = 'D_PBL'; val = 'PBL'; }, { key = 'D_CAPE'; val = 'CAPE'; } + { key = 'D_MLCAPE'; val = 'MLCAPE'; } ]; diff --git a/met/scripts/config/EnsembleStatConfig b/met/scripts/config/EnsembleStatConfig index 91d4cab5d0..64367110ea 100644 --- a/met/scripts/config/EnsembleStatConfig +++ b/met/scripts/config/EnsembleStatConfig @@ -249,6 +249,10 @@ interp = { { method = NEAREST; width = 1; + }, + { + method = HIRA; + width = 2; } ]; } diff --git a/met/scripts/config/PB2NCConfig_G212 b/met/scripts/config/PB2NCConfig_G212 index d5fbe9c141..0edaddf6e6 100644 --- a/met/scripts/config/PB2NCConfig_G212 +++ b/met/scripts/config/PB2NCConfig_G212 @@ -84,7 +84,7 @@ obs_bufr_var = [ "QOB", "TOB", "ZOB", "UOB", "VOB", // Mapping of BUFR variable name to GRIB name. The default map is defined at // obs_prepbufr_map. This replaces/expends the default map. // -//obs_bufr_var = []; +obs_bufr_map = []; //////////////////////////////////////////////////////////////////////////////// diff --git a/met/scripts/config/STATAnalysisConfig b/met/scripts/config/STATAnalysisConfig index 5f81208f5a..712813c57e 100644 --- a/met/scripts/config/STATAnalysisConfig +++ b/met/scripts/config/STATAnalysisConfig @@ -84,7 +84,7 @@ jobs = [ "-job aggregate -line_type GRAD -vx_mask DTC165 -vx_mask DTC166 -by FCST_VAR -dump_row ${TEST_OUT_DIR}/stat_analysis/job_aggregate_GRAD.stat", "-job aggregate -line_type ISC -fcst_thresh >0.0 -vx_mask TILE_TOT -fcst_var APCP_12 -dump_row ${TEST_OUT_DIR}/stat_analysis/job_aggregate_ISC.stat", "-job aggregate -line_type RHIST -obtype MC_PCP -vx_mask HUC4_1605 -vx_mask HUC4_1803 -vx_mask HUC4_1804 -vx_mask HUC4_1805 -vx_mask HUC4_1806 -dump_row ${TEST_OUT_DIR}/stat_analysis/job_aggregate_RHIST.stat", - "-job aggregate_stat -line_type ORANK -out_line_type RHIST -obtype ADPSFC -vx_mask HUC4_1605 -vx_mask HUC4_1803 -vx_mask HUC4_1804 -vx_mask HUC4_1805 -vx_mask HUC4_1806 -dump_row ${TEST_OUT_DIR}/stat_analysis/job_aggregate_stat_ORANK_RHIST.stat" + "-job aggregate_stat -line_type ORANK -out_line_type RHIST -obtype ADPSFC -vx_mask HUC4_1605,HUC4_1803,HUC4_1804,HUC4_1805,HUC4_1806 -by INTERP_MTHD,INTERP_PNTS -dump_row ${TEST_OUT_DIR}/stat_analysis/job_aggregate_stat_ORANK_RHIST.stat" ]; //////////////////////////////////////////////////////////////////////////////// diff --git a/met/src/basic/vx_config/config_constants.h b/met/src/basic/vx_config/config_constants.h index 4649c555a4..f3c5550f02 100644 --- a/met/src/basic/vx_config/config_constants.h +++ b/met/src/basic/vx_config/config_constants.h @@ -571,7 +571,7 @@ static const char conf_key_message_type_group_map[] = "message_type_group_map"; static const char conf_key_obs_bufr_map[] = "obs_bufr_map"; static const char conf_key_obs_bufr_var[] = "obs_bufr_var"; static const char conf_key_obs_name_map[] = "obs_name_map"; -static const char conf_key_obs_prefbufr_map[] = "obs_prefbufr_map"; +static const char conf_key_obs_prepbufr_map[] = "obs_prepbufr_map"; static const char conf_key_key[] = "key"; static const char conf_key_val[] = "val"; static const char conf_key_boot_interval[] = "boot.interval"; diff --git a/met/src/basic/vx_config/config_util.cc b/met/src/basic/vx_config/config_util.cc index bfdeef07b2..edeabb649a 100644 --- a/met/src/basic/vx_config/config_util.cc +++ b/met/src/basic/vx_config/config_util.cc @@ -24,6 +24,7 @@ using namespace std; /////////////////////////////////////////////////////////////////////////////// static const double default_vld_thresh = 1.0; +static const char conf_key_prepbufr_map_bad[] = "obs_prefbufr_map"; // for backward compatibility /////////////////////////////////////////////////////////////////////////////// @@ -140,9 +141,10 @@ RegridInfo::RegridInfo() { void RegridInfo::validate() { // Check for unsupported regridding options - if(method == InterpMthd_Best || + if(method == InterpMthd_Best || method == InterpMthd_Geog_Match || - method == InterpMthd_Gaussian) { + method == InterpMthd_Gaussian || + method == InterpMthd_HiRA) { mlog << Error << "\nRegridInfo::validate() -> " << "\"" << interpmthd_to_string(method) << "\" not valid for regridding, only interpolating.\n\n"; @@ -1088,14 +1090,6 @@ map parse_conf_metadata_map(Dictionary *dict) { /////////////////////////////////////////////////////////////////////////////// -map parse_conf_obs_bufr_map(Dictionary *dict) { - map m = parse_conf_key_value_map(dict, conf_key_obs_prefbufr_map); - parse_add_conf_key_value_map(dict, conf_key_obs_bufr_map, &m); - return m; -} - -/////////////////////////////////////////////////////////////////////////////// - map parse_conf_obs_name_map(Dictionary *dict) { const char *method_name = "parse_conf_obs_name_map() -> "; return parse_conf_key_value_map(dict, conf_key_obs_name_map); @@ -2279,6 +2273,7 @@ InterpMthd int_to_interpmthd(int i) { else if(i == conf_const.lookup_int(interpmthd_gaussian_str)) m = InterpMthd_Gaussian; else if(i == conf_const.lookup_int(interpmthd_maxgauss_str)) m = InterpMthd_MaxGauss; else if(i == conf_const.lookup_int(interpmthd_geog_match_str)) m = InterpMthd_Geog_Match; + else if(i == conf_const.lookup_int(interpmthd_hira_str)) m = InterpMthd_HiRA; else { mlog << Error << "\nconf_int_to_interpmthd() -> " << "Unexpected value of " << i diff --git a/met/src/basic/vx_config/config_util.h b/met/src/basic/vx_config/config_util.h index f9750e0212..5afcc974a1 100644 --- a/met/src/basic/vx_config/config_util.h +++ b/met/src/basic/vx_config/config_util.h @@ -25,6 +25,10 @@ //////////////////////////////////////////////////////////////////////// +static const char conf_key_old_prepbufr_map[] = "obs_prefbufr_map"; // for backward compatibility + +//////////////////////////////////////////////////////////////////////// + extern ConcatString parse_conf_version(Dictionary *dict); extern ConcatString parse_conf_string(Dictionary *dict, const char *, bool check_empty = true); extern GrdFileType parse_conf_file_type(Dictionary *dict); @@ -48,13 +52,13 @@ extern ClimoCDFInfo parse_conf_climo_cdf(Dictionary *dict); extern TimeSummaryInfo parse_conf_time_summary(Dictionary *dict); extern map parse_conf_key_value_map( Dictionary *dict, const char *conf_key_map_name, const char *caller=0); +extern void parse_add_conf_key_value_map( + Dictionary *dict, const char *conf_key_map_name, map *m); extern map parse_conf_message_type_map(Dictionary *dict); extern map parse_conf_message_type_group_map(Dictionary *dict); extern map parse_conf_metadata_map(Dictionary *dict); -extern map - parse_conf_obs_bufr_map(Dictionary *dict); extern map parse_conf_obs_name_map(Dictionary *dict); extern BootInfo parse_conf_boot(Dictionary *dict); diff --git a/met/src/basic/vx_util/interp_mthd.cc b/met/src/basic/vx_util/interp_mthd.cc index dcd9ff5ac4..d84bf9c4ab 100644 --- a/met/src/basic/vx_util/interp_mthd.cc +++ b/met/src/basic/vx_util/interp_mthd.cc @@ -44,6 +44,7 @@ ConcatString interpmthd_to_string(const InterpMthd m) { case(InterpMthd_Gaussian): out = interpmthd_gaussian_str; break; case(InterpMthd_MaxGauss): out = interpmthd_maxgauss_str; break; case(InterpMthd_Geog_Match): out = interpmthd_geog_match_str; break; + case(InterpMthd_HiRA): out = interpmthd_hira_str; break; case(InterpMthd_None): default: out = interpmthd_none_str; break; @@ -77,6 +78,7 @@ InterpMthd string_to_interpmthd(const char *mthd_str) { else if(strcmp(mthd_str, interpmthd_gaussian_str ) == 0) m = InterpMthd_Gaussian; else if(strcmp(mthd_str, interpmthd_maxgauss_str ) == 0) m = InterpMthd_MaxGauss; else if(strcmp(mthd_str, interpmthd_geog_match_str) == 0) m = InterpMthd_Geog_Match; + else if(strcmp(mthd_str, interpmthd_hira_str) == 0) m = InterpMthd_HiRA; else m = InterpMthd_None; return(m); diff --git a/met/src/basic/vx_util/interp_mthd.h b/met/src/basic/vx_util/interp_mthd.h index e65be8f7ef..39c4b23e87 100644 --- a/met/src/basic/vx_util/interp_mthd.h +++ b/met/src/basic/vx_util/interp_mthd.h @@ -42,7 +42,8 @@ enum InterpMthd { InterpMthd_Lower_Left, InterpMthd_Gaussian, InterpMthd_MaxGauss, - InterpMthd_Geog_Match + InterpMthd_Geog_Match, + InterpMthd_HiRA }; // @@ -69,6 +70,7 @@ static const char interpmthd_lower_left_str[] = "LOWER_LEFT"; static const char interpmthd_gaussian_str[] = "GAUSSIAN"; static const char interpmthd_maxgauss_str[] = "MAXGAUSS"; static const char interpmthd_geog_match_str[] = "GEOG_MATCH"; +static const char interpmthd_hira_str[] = "HIRA"; /////////////////////////////////////////////////////////////////////////////// diff --git a/met/src/basic/vx_util/interp_util.cc b/met/src/basic/vx_util/interp_util.cc index 631b88d265..bdaf7a573d 100644 --- a/met/src/basic/vx_util/interp_util.cc +++ b/met/src/basic/vx_util/interp_util.cc @@ -1020,6 +1020,7 @@ double compute_sfc_interp(const DataPlane &dp, mlog << Error << "\ncompute_sfc_interp() -> " << "unsupported interpolation method encountered: " << interpmthd_to_string(mthd) << "(" << mthd << ")\n\n"; + exit(1); } delete gt; diff --git a/met/src/libcode/vx_data2d/var_info.h b/met/src/libcode/vx_data2d/var_info.h index 382f25aded..3661b07ee3 100644 --- a/met/src/libcode/vx_data2d/var_info.h +++ b/met/src/libcode/vx_data2d/var_info.h @@ -258,9 +258,9 @@ inline int VarInfo::accum_attr() const { return(SetAttrAccum); } // struct InputInfo { - VarInfo * var_info; // Variable information to read - int file_index; // Index in file_list of file to read - StringArray * file_list; // Array of files (unallocated) + VarInfo * var_info; // Variable information to read + int file_index; // Index in file_list of file to read + StringArray * file_list; // Array of files (unallocated) }; //////////////////////////////////////////////////////////////////////// @@ -270,33 +270,37 @@ struct InputInfo { // class EnsVarInfo { -private: - vector inputs; // Vector of InputInfo - VarInfo * ctrl_info; // Field info for control member -public: - EnsVarInfo(); - ~EnsVarInfo(); - EnsVarInfo(const EnsVarInfo &); + private: + vector inputs; // Vector of InputInfo + VarInfo * ctrl_info; // Field info for control member - void clear(); - void assign(const EnsVarInfo &); + public: + EnsVarInfo(); + ~EnsVarInfo(); + EnsVarInfo(const EnsVarInfo &); + + void clear(); + void assign(const EnsVarInfo &); - void add_input(InputInfo); - int inputs_n(); + void add_input(InputInfo); + int inputs_n(); - void set_ctrl(VarInfo *); - VarInfo * get_ctrl(int); + void set_ctrl(VarInfo *); + VarInfo * get_ctrl(int); - // Get VarInfo from first InputInfo if requested VarInfo is NULL - VarInfo * get_var_info(int index=0); - ConcatString get_file(int index=0); - int get_file_index(int index=0); + // Get VarInfo from first InputInfo if requested VarInfo is NULL + VarInfo * get_var_info(int index=0); + ConcatString get_file(int index=0); + int get_file_index(int index=0); + + ConcatString nc_var_str; // Ensemble variable name strings + ThreshArray cat_ta; // Ensemble categorical thresholds + ConcatString raw_magic_str; // Magic string w/o var substitution - ConcatString nc_var_str; // Ensemble variable name strings - ThreshArray cat_ta; // Ensemble categorical thresholds - ConcatString raw_magic_str; // Magic string w/o var substitution }; +//////////////////////////////////////////////////////////////////////// + ConcatString raw_magic_str(Dictionary i_edict, GrdFileType file_type); /////////////////////////////////////////////////////////////////////////////// diff --git a/met/src/libcode/vx_data2d_nc_met/met_file.cc b/met/src/libcode/vx_data2d_nc_met/met_file.cc index dd8f599c48..14a41b64d1 100644 --- a/met/src/libcode/vx_data2d_nc_met/met_file.cc +++ b/met/src/libcode/vx_data2d_nc_met/met_file.cc @@ -330,8 +330,6 @@ void MetNcFile::dump(ostream & out, int depth) const { int j, k; -int month, day, year, hour, minute, second; -char junk[256]; Indent prefix(depth); Indent p2(depth + 1); Indent p3(depth + 2); @@ -362,12 +360,6 @@ out << prefix << "Ydim = " << (Ydim ? GET_NC_NAME_P(Ydim) : "(nul)") << "\n"; out << prefix << "\n"; -snprintf(junk, sizeof(junk), "%s %d, %d %2d:%02d:%02d", short_month_name[month], day, year, hour, minute, second); - -out << junk << "\n"; - -out << prefix << "\n"; - out << prefix << "Nvars = " << Nvars << "\n"; for (j=0; j " @@ -1276,7 +1280,7 @@ void VxPairDataEnsemble::set_interp(int i_interp, //////////////////////////////////////////////////////////////////////// void VxPairDataEnsemble::set_interp(int i_interp, InterpMthd mthd, - int width, GridTemplateFactory::GridTemplates shape) { + int width, GridTemplateFactory::GridTemplates shape) { for(int i=0; isize()); + } + else { + pd[i][j][k].set_ens_size(n); + } } } } @@ -1624,9 +1639,10 @@ void VxPairDataEnsemble::add_point_obs(float *hdr_arr, int *hdr_typ_arr, //////////////////////////////////////////////////////////////////////// void VxPairDataEnsemble::add_ens(int member, bool mn, Grid &gr) { - int i, j, k, l; - int f_lvl_blw, f_lvl_abv; + int i, j, k, l, m; + int f_lvl_blw, f_lvl_abv, i_mem; double to_lvl, fcst_v; + NumArray fcst_na; // Set flag for specific humidity bool spfh_flag = fcst_info->get_var_info()->is_specific_humidity() && @@ -1637,9 +1653,24 @@ void VxPairDataEnsemble::add_ens(int member, bool mn, Grid &gr) { for(j=0; j " + << "the \"" << interpmthd_to_string(pd[0][0][k].interp_mthd) + << "\" interpolation method only applies when verifying a " + << "single level, not " << fcst_dpa.n_planes() + << " levels.\n\n"; + continue; + } + // Process each of the observations for(l=0; lget_var_info()->level().type() == LevelType_Pres ? pd[i][j][k].lvl_na[l] : pd[i][j][k].elv_na[l]); @@ -1655,46 +1686,83 @@ void VxPairDataEnsemble::add_ens(int member, bool mn, Grid &gr) { find_vert_lvl(fcst_dpa, to_lvl, f_lvl_blw, f_lvl_abv); } - // Compute the interpolated ensemble value - fcst_v = compute_interp(fcst_dpa, - pd[i][j][k].x_na[l], - pd[i][j][k].y_na[l], - pd[i][j][k].o_na[l], - pd[i][j][k].cmn_na[l], - pd[i][j][k].csd_na[l], - pd[0][0][k].interp_mthd, - pd[0][0][k].interp_wdth, - pd[0][0][k].interp_shape, - gr.wrap_lon(), - interp_thresh, spfh_flag, - fcst_info->get_var_info()->level().type(), - to_lvl, f_lvl_blw, f_lvl_abv); - - // Store the ensemble mean - if(mn) { - pd[i][j][k].mn_na.add(fcst_v); + // Extract the HiRA neighborhood of values + if(pd[0][0][k].interp_mthd == InterpMthd_HiRA) { + + // For HiRA, set the ensemble mean to bad data + if(mn) { + fcst_na.erase(); + fcst_na.add(bad_data_double); + } + // Otherwise, retrieve all the neighborhood values + // using a valid threshold of 0 + else { + get_interp_points(fcst_dpa, + pd[i][j][k].x_na[l], + pd[i][j][k].y_na[l], + pd[0][0][k].interp_mthd, + pd[0][0][k].interp_wdth, + pd[0][0][k].interp_shape, + gr.wrap_lon(), + 0, spfh_flag, + fcst_info->get_var_info()->level().type(), + to_lvl, f_lvl_blw, f_lvl_abv, + fcst_na); + } } - // Store the ensemble member values + // Otherwise, get a single interpolated ensemble value else { + fcst_na.add(compute_interp(fcst_dpa, + pd[i][j][k].x_na[l], + pd[i][j][k].y_na[l], + pd[i][j][k].o_na[l], + pd[i][j][k].cmn_na[l], + pd[i][j][k].csd_na[l], + pd[0][0][k].interp_mthd, + pd[0][0][k].interp_wdth, + pd[0][0][k].interp_shape, + gr.wrap_lon(), + interp_thresh, spfh_flag, + fcst_info->get_var_info()->level().type(), + to_lvl, f_lvl_blw, f_lvl_abv)); + } - // Track unperturbed ensemble variance sums - // Exclude the control member from the variance - if(member != pd[i][j][k].ctrl_index) { - pd[i][j][k].add_ens_var_sums(l, fcst_v); - } + // Store the single ensemble value or HiRA neighborhood + for(m=0; mflag) { - fcst_v = add_obs_error_inc( - obs_error_info->rng_ptr, FieldType_Fcst, - pd[i][j][k].obs_error_entry[l], - pd[i][j][k].o_na[l], fcst_v); + // Store the ensemble mean + if(mn) { + pd[i][j][k].mn_na.add(fcst_na[m]); + } + // Store the ensemble member values + else { + + // Track unperturbed ensemble variance sums + // Exclude the control member from the variance + if(member != pd[i][j][k].ctrl_index) { + pd[i][j][k].add_ens_var_sums(l, fcst_na[m]); + } + + // Apply observation error perturbation, if requested + if(obs_error_info->flag) { + fcst_v = add_obs_error_inc( + obs_error_info->rng_ptr, FieldType_Fcst, + pd[i][j][k].obs_error_entry[l], + pd[i][j][k].o_na[l], fcst_na[m]); + } + else { + fcst_v = fcst_na[m]; + } + + // Determine index of ensemble member + i_mem = member * fcst_na.n() + m; + + // Store perturbed ensemble member value + pd[i][j][k].add_ens(i_mem, fcst_v); } - // Store perturbed ensemble member value - pd[i][j][k].add_ens(member, fcst_v); - } - } + } // end for m - fcst_na + } // end for l - n_obs } // end for k - n_interp } // end for j - n_mask } // end for i - n_msg_typ diff --git a/met/src/tools/core/ensemble_stat/ensemble_stat.cc b/met/src/tools/core/ensemble_stat/ensemble_stat.cc index 078c547e8a..ce0784e84e 100644 --- a/met/src/tools/core/ensemble_stat/ensemble_stat.cc +++ b/met/src/tools/core/ensemble_stat/ensemble_stat.cc @@ -64,6 +64,7 @@ // 032 10/07/21 Halley Gotway MET #1905 Add -ctrl option. // 033 11/15/21 Halley Gotway MET #1968 Ensemble -ctrl error check. // 034 01/14/21 McCabe MET #1695 All members in one file. +// 035 02/15/22 Halley Gotway MET #1583 Add HiRA option. // //////////////////////////////////////////////////////////////////////// @@ -617,15 +618,20 @@ void process_n_vld() { n_vx_vld.clear(); // Loop through the ensemble fields to be processed - var_it = conf_info.ens_input.begin(); - n_ens_inputs = (*var_it)->inputs_n(); + if(conf_info.ens_input.size() > 0) { + var_it = conf_info.ens_input.begin(); + n_ens_inputs = (*var_it)->inputs_n(); + } + else { + n_ens_inputs = 0; + } for(i_var=0; var_it != conf_info.ens_input.end(); var_it++, i_var++) { // Loop through the ensemble inputs for(i_ens=n_vld=0; i_ens < n_ens_inputs; i_ens++) { - // get file and VarInfo to process + // Get file and VarInfo to process ens_file = (*var_it)->get_file(i_ens); var_info = (*var_it)->get_var_info(i_ens); @@ -1276,7 +1282,7 @@ int process_point_ens(int i_ens, int &n_miss) { info = conf_info.vx_opt[i].vx_pd.fcst_info->get_var_info(i_ens); } - // if not processing mean, get file based on current vx and ensemble index + // If not processing mean, get file based on current vx and ensemble index if(!is_ens_mean) ens_file = conf_info.vx_opt[i].vx_pd.fcst_info->get_file(i_ens); // Read the gridded data from the input forecast file @@ -1381,9 +1387,11 @@ void process_point_scores() { for(l=0; l " - << mthd_str << " smoothing option not supported for " - << "gridded observations.\n\n"; + << mthd_str << " option not supported for " + << "smoothing gridded observations.\n\n"; continue; } @@ -2306,7 +2314,7 @@ void setup_nc_file(unixtime valid_ut, const char *suffix) { //////////////////////////////////////////////////////////////////////// void setup_txt_files() { - int i, n, n_phist_bin, n_vld, max_col; + int i, n, n_phist_bin, max_n_ens, max_col; ConcatString tmp_str; // Check to see if the text files have already been set up @@ -2327,16 +2335,19 @@ void setup_txt_files() { n_phist_bin = (n > n_phist_bin ? n : n_phist_bin); } - // Store the maximum number of valid verification fields - n_vld = n_vx_vld.max(); + // Get the maximum number of ensemble members, including HiRA + max_n_ens = n_vx_vld.max(); + if(conf_info.get_max_hira_size() > 0) { + max_n_ens *= conf_info.get_max_hira_size(); + } // Get the maximum number of data columns - max_col = max(get_n_orank_columns(n_vld+1), - get_n_rhist_columns(n_vld)); + max_col = max(get_n_orank_columns(max_n_ens+1), + get_n_rhist_columns(max_n_ens)); max_col = max(max_col, get_n_phist_columns(n_phist_bin)); max_col = max(max_col, n_ecnt_columns); max_col = max(max_col, n_ssvar_columns); - max_col = max(max_col, get_n_relp_columns(n_vld)); + max_col = max(max_col, get_n_relp_columns(max_n_ens)); max_col += n_header_columns; // Initialize file stream @@ -2387,7 +2398,7 @@ void setup_txt_files() { switch(i) { case(i_rhist): - max_col = get_n_rhist_columns(n_vld+1) + n_header_columns + 1; + max_col = get_n_rhist_columns(max_n_ens+1) + n_header_columns + 1; break; case(i_phist): @@ -2395,11 +2406,11 @@ void setup_txt_files() { break; case(i_relp): - max_col = get_n_relp_columns(n_vld) + n_header_columns + 1; + max_col = get_n_relp_columns(max_n_ens) + n_header_columns + 1; break; case(i_orank): - max_col = get_n_orank_columns(n_vld) + n_header_columns + 1; + max_col = get_n_orank_columns(max_n_ens) + n_header_columns + 1; break; default: @@ -2415,7 +2426,7 @@ void setup_txt_files() { switch(i) { case(i_rhist): - write_rhist_header_row(1, n_vld+1, txt_at[i], 0, 0); + write_rhist_header_row(1, max_n_ens+1, txt_at[i], 0, 0); break; case(i_phist): @@ -2423,11 +2434,11 @@ void setup_txt_files() { break; case(i_relp): - write_relp_header_row(1, n_vld, txt_at[i], 0, 0); + write_relp_header_row(1, max_n_ens, txt_at[i], 0, 0); break; case(i_orank): - write_orank_header_row(1, n_vld, txt_at[i], 0, 0); + write_orank_header_row(1, max_n_ens, txt_at[i], 0, 0); break; default: diff --git a/met/src/tools/core/ensemble_stat/ensemble_stat_conf_info.cc b/met/src/tools/core/ensemble_stat/ensemble_stat_conf_info.cc index e9aa6a2765..252d7eaeb9 100644 --- a/met/src/tools/core/ensemble_stat/ensemble_stat_conf_info.cc +++ b/met/src/tools/core/ensemble_stat/ensemble_stat_conf_info.cc @@ -98,10 +98,11 @@ void EnsembleStatConfInfo::clear() { ens_input.clear(); // Reset counts - n_ens_var = 0; - max_n_thresh = 0; - n_nbrhd = 0; - n_vx = 0; + n_ens_var = 0; + n_nbrhd = 0; + n_vx = 0; + max_n_thresh = 0; + max_hira_size = 0; return; } @@ -132,7 +133,7 @@ void EnsembleStatConfInfo::process_config(GrdFileType etype, StringArray * ens_files, StringArray * fcst_files, bool use_ctrl) { - int i,j, n_ens_files; + int i, j, n_ens_files; VarInfoFactory info_factory; mapoutput_map; Dictionary *edict = (Dictionary *) 0; @@ -230,14 +231,11 @@ void EnsembleStatConfInfo::process_config(GrdFileType etype, } // If no ensemble member IDs were provided, add an empty string - if(ens_member_ids.n() == 0) { - ens_member_ids.add(""); - } + if(ens_member_ids.n() == 0) ens_member_ids.add(""); // Conf: ens.field edict = conf.lookup_array(conf_key_ens_field); - // Determine the number of ensemble fields to be processed n_ens_var = parse_conf_n_vx(edict); @@ -251,7 +249,7 @@ void EnsembleStatConfInfo::process_config(GrdFileType etype, } // Parse the ensemble field information - for(i=0,max_n_thresh=0; iset_dict(i_edict); - // Dump the contents of the current VarInfo if(mlog.verbosity_level() >= 5) { mlog << Debug(5) @@ -296,7 +293,6 @@ void EnsembleStatConfInfo::process_config(GrdFileType etype, input_info.file_list = ens_files; ens_info->add_input(input_info); } // end for k - } // end for j // Get field info for control member if set @@ -315,7 +311,7 @@ void EnsembleStatConfInfo::process_config(GrdFileType etype, } // Conf: ens_nc_var_str - ens_info->nc_var_str =parse_conf_string(&i_edict, conf_key_nc_var_str, false); + ens_info->nc_var_str = parse_conf_string(&i_edict, conf_key_nc_var_str, false); // Conf: ens_nc_pairs // Only parse thresholds if probabilities are requested @@ -332,8 +328,8 @@ void EnsembleStatConfInfo::process_config(GrdFileType etype, } // Keep track of the maximum number of thresholds - if(ens_info->cat_ta.n_elements() > max_n_thresh) { - max_n_thresh = ens_info->cat_ta.n_elements(); + if(ens_info->cat_ta.n() > max_n_thresh) { + max_n_thresh = ens_info->cat_ta.n(); } } @@ -366,11 +362,11 @@ void EnsembleStatConfInfo::process_config(GrdFileType etype, } // Conf: nbrhd_prob - nbrhd_prob = parse_conf_nbrhd(edict, conf_key_nbrhd_prob); + nbrhd_prob = parse_conf_nbrhd(&conf, conf_key_nbrhd_prob); n_nbrhd = nbrhd_prob.width.n(); // Conf: nmep_smooth - nmep_smooth = parse_conf_interp(edict, conf_key_nmep_smooth); + nmep_smooth = parse_conf_interp(&conf, conf_key_nmep_smooth); // Loop through the neighborhood probability smoothing options for(i=0; isize()); + } + } } // Summarize output flags across all verification tasks @@ -580,7 +587,7 @@ void EnsembleStatConfInfo::process_masks(const Grid &grid) { vx_opt[i].mask_name_area.clear(); // Parse the masking grids - for(j=0; j " << "At least one grid, polyline or station ID masking " << "region must be provided for verification task number " @@ -740,15 +747,18 @@ void EnsembleStatVxOpt::clear() { vx_pd.clear(); var_str.clear(); beg_ds = end_ds = bad_data_int; + mask_grid.clear(); mask_poly.clear(); mask_sid.clear(); mask_llpnt.clear(); mask_name.clear(); mask_name_area.clear(); + msg_typ.clear(); othr_ta.clear(); cdf_info.clear(); + ci_alpha.clear(); interp_info.clear(); @@ -992,9 +1002,9 @@ void EnsembleStatVxOpt::process_config(GrdFileType ftype, Dictionary &fdict, //////////////////////////////////////////////////////////////////////// void EnsembleStatVxOpt::set_vx_pd(EnsembleStatConfInfo *conf_info, int ctrl_index) { - int i, n; - int n_msg_typ = msg_typ.n_elements(); - int n_mask = mask_name.n_elements(); + int i, j, n; + int n_msg_typ = msg_typ.n(); + int n_mask = mask_name.n(); int n_interp = interp_info.n_interp; StringArray sa; @@ -1041,41 +1051,41 @@ void EnsembleStatVxOpt::set_vx_pd(EnsembleStatConfInfo *conf_info, int ctrl_inde for(i=0; imsg_typ_group_map[msg_typ[i]]; - if(sa.n_elements() == 0) sa.add(msg_typ[i]); + if(sa.n() == 0) sa.add(msg_typ[i]); vx_pd.set_msg_typ_vals(i, sa); } // Define the masking information: grid, poly, sid // Define the grid masks - for(i=0; imask_area_map[mask_name[n]])); } // Define the poly masks - for(i=0; imask_area_map[mask_name[n]])); } // Define the station ID masks - for(i=0; imask_sid_map[mask_name[n]])); } // Define the Lat/Lon point masks for(i=0; i<(int) mask_llpnt.size(); i++) { - n = i + mask_grid.n_elements() + mask_poly.n_elements() + mask_sid.n_elements(); + n = i + mask_grid.n() + mask_poly.n() + mask_sid.n(); vx_pd.set_mask_llpnt(n, mask_name[n].c_str(), &mask_llpnt[i]); } // Define the interpolation methods - for(i=0; i ens_input; // Vector of EnsVarInfo pointers (allocated) StringArray ens_member_ids; // Array of ensemble member ID strings ConcatString control_id; // Control ID - NbrhdInfo nbrhd_prob; // Neighborhood probability definition - int n_nbrhd; // Number of neighborhood sizes - InterpInfo nmep_smooth; // Neighborhood maximum smoothing information + NbrhdInfo nbrhd_prob; // Neighborhood probability definition + int n_nbrhd; // Number of neighborhood sizes + InterpInfo nmep_smooth; // Neighborhood maximum smoothing information - EnsembleStatVxOpt * vx_opt; // Array of vx task options [n_vx] (allocated) + EnsembleStatVxOpt * vx_opt; // Array of vx task options [n_vx] (allocated) - double vld_ens_thresh; // Required ratio of valid input files - double vld_data_thresh; // Required ratio of valid data for each point + double vld_ens_thresh; // Required ratio of valid input files + double vld_data_thresh; // Required ratio of valid data for each point // Message type groups that should be processed together map msg_typ_group_map; @@ -245,25 +246,30 @@ class EnsembleStatConfInfo { void set_vx_pd (const IntArray &, int); // Dump out the counts - int get_n_ens_var() const; - int get_max_n_thresh() const; - int get_n_nbrhd() const; - int get_n_vx() const; + int get_n_ens_var() const; + int get_n_nbrhd() const; + int get_n_vx() const; // Compute the maximum number of output lines possible based // on the contents of the configuration file int n_txt_row(int i) const; int n_stat_row() const; + + // Maximum across all verification tasks + int get_max_n_thresh() const; + int get_max_hira_size() const; + int get_compression_level(); }; //////////////////////////////////////////////////////////////////////// -inline int EnsembleStatConfInfo::get_n_ens_var() const { return(n_ens_var); } -inline int EnsembleStatConfInfo::get_max_n_thresh() const { return(max_n_thresh); } -inline int EnsembleStatConfInfo::get_n_nbrhd() const { return(n_nbrhd); } -inline int EnsembleStatConfInfo::get_n_vx() const { return(n_vx); } -inline int EnsembleStatConfInfo::get_compression_level() { return(conf.nc_compression()); } +inline int EnsembleStatConfInfo::get_n_ens_var() const { return(n_ens_var); } +inline int EnsembleStatConfInfo::get_n_nbrhd() const { return(n_nbrhd); } +inline int EnsembleStatConfInfo::get_n_vx() const { return(n_vx); } +inline int EnsembleStatConfInfo::get_max_n_thresh() const { return(max_n_thresh); } +inline int EnsembleStatConfInfo::get_max_hira_size() const { return(max_hira_size); } +inline int EnsembleStatConfInfo::get_compression_level() { return(conf.nc_compression()); } //////////////////////////////////////////////////////////////////////// diff --git a/met/src/tools/core/point_stat/point_stat.cc b/met/src/tools/core/point_stat/point_stat.cc index edc3c02f1a..73dfa2af24 100644 --- a/met/src/tools/core/point_stat/point_stat.cc +++ b/met/src/tools/core/point_stat/point_stat.cc @@ -1221,7 +1221,7 @@ void process_scores() { // Appy HiRA verification and write ensemble output do_hira_ens(i, pd_ptr); - } // end HiRA for probabilities + } // end HiRA for ensembles // Apply HiRA probabilistic verification logic if(!conf_info.vx_opt[i].vx_pd.fcst_info->is_prob() && diff --git a/met/src/tools/core/point_stat/point_stat_conf_info.cc b/met/src/tools/core/point_stat/point_stat_conf_info.cc index 2d639a20d6..49f1de8a9c 100644 --- a/met/src/tools/core/point_stat/point_stat_conf_info.cc +++ b/met/src/tools/core/point_stat/point_stat_conf_info.cc @@ -1189,8 +1189,7 @@ int PointStatVxOpt::n_txt_row(int i_txt_row) const { n = (!prob_flag ? 0 : n_pd * get_n_oprob_thresh() * n_bin); // Number of HiRA PCT, PJC, or PRC lines = - // Message Types * Masks * Interpolations * Thresholds * - // HiRA widths + // Message Types * Masks * HiRA widths * Thresholds if(hira_info.flag) { n += (prob_flag ? 0 : n_pd * get_n_cat_thresh() * hira_info.width.n()); @@ -1206,8 +1205,8 @@ int PointStatVxOpt::n_txt_row(int i_txt_row) const { get_n_oprob_thresh() * get_n_ci_alpha() * n_bin); // Number of HiRA PSTD lines = - // Message Types * Masks * Interpolations * Thresholds * - // HiRA widths * Alphas + // Message Types * Masks * HiRA widths * Thresholds * + // Alphas if(hira_info.flag) { n += (prob_flag ? 0 : n_pd * get_n_cat_thresh() * hira_info.width.n() * @@ -1218,8 +1217,8 @@ int PointStatVxOpt::n_txt_row(int i_txt_row) const { case(i_ecnt): case(i_rps): - // Number of HiRA ECNT, ORANK and RPS lines = - // Message Types * Masks * Interpolations * HiRA widths * + // Number of HiRA ECNT and RPS lines = + // Message Types * Masks * HiRA widths * // Alphas if(hira_info.flag) { n = n_pd * hira_info.width.n() * get_n_ci_alpha(); @@ -1231,7 +1230,7 @@ int PointStatVxOpt::n_txt_row(int i_txt_row) const { break; case(i_orank): - // Number HiRA ORANK lines possible = + // Number of HiRA ORANK lines possible = // Number of pairs * Categorical Thresholds * // HiRA widths if(hira_info.flag) { diff --git a/met/src/tools/core/series_analysis/series_analysis.cc b/met/src/tools/core/series_analysis/series_analysis.cc index 6a33eefd73..72a897fa02 100644 --- a/met/src/tools/core/series_analysis/series_analysis.cc +++ b/met/src/tools/core/series_analysis/series_analysis.cc @@ -670,9 +670,10 @@ void process_scores() { // Number of points skipped due to valid data threshold int n_skip_zero = 0; int n_skip_pos = 0; - + // Loop over the data reads for(i_read=0; i_read: XNEW(1)',XOLD(1),XNEW(1) -C print*,' DEBUG HS calcape XOLD(2)->: XNEW(2)',XOLD(2),XNEW(2) -CC print*,' DEBUG HS calcape XNEW: ',XNEW -CC print*,' DEBUG HS calcape YNEW: ',YNEW -C print*,' DEBUG HS calcape YOLD(1)->: YNEW(1)',YOLD(1),YNEW(1) -C print*,' DEBUG HS calcape YOLD(2)->: YNEW(2)',YOLD(2),YNEW(2) -C print*,' DEBUG HS calcape YOLD(3)->: YNEW(3)',YOLD(3),YNEW(3) -C print*,' DEBUG HS calcape YOLD(-2)->: YNEW(-2)', C 2 YOLD(NOLD-1),YNEW(NNEW-1) -C print*,' DEBUG HS calcape YOLD(-1)->: YNEW(-1)', C 2 YOLD(NOLD),YNEW(NNEW) RETURN END diff --git a/met/src/tools/other/pb2nc/pb2nc.cc b/met/src/tools/other/pb2nc/pb2nc.cc index 7df4d754eb..f08fca1038 100644 --- a/met/src/tools/other/pb2nc/pb2nc.cc +++ b/met/src/tools/other/pb2nc/pb2nc.cc @@ -260,8 +260,10 @@ static const char *airnow_aux_vars = "TPHR QCIND"; // Pick the latter one if exists multiuple variables static const char *bufr_avail_sid_names = "SID SAID RPID"; static const char *bufr_avail_latlon_names = "XOB CLON CLONH YOB CLAT CLATH"; +static const char *derived_mlcape = "D_MLCAPE"; static const char *derived_cape = "D_CAPE"; static const char *derived_pbl = "D_PBL"; +static const float MLCAPE_INTERVAL = 3000.; static double bufr_obs[mxr8lv][mxr8pm]; static double bufr_obs_extra[mxr8lv][mxr8pm]; @@ -511,6 +513,7 @@ void initialize() { prepbufr_derive_vars.add("D_MIXR"); prepbufr_derive_vars.add("D_PRMSL"); prepbufr_derive_vars.add("D_CAPE"); + prepbufr_derive_vars.add("D_MLCAPE"); prepbufr_derive_vars.add("D_PBL"); for (int idx=0; idx<(sizeof(hdr) / sizeof(hdr[0])); idx++) { @@ -982,10 +985,13 @@ void process_pbfile(int i_pb) { int itype = 1; // itype 1: where a parcel is lifted from the ground // itype 2: Where the "best cape" in a number of parcels int cape_code = -1; + int mlcape_code = -1; float p1d,t1d,q1d; int IMM, JMM; - int cape_level=0, cape_count=0, cape_cnt_too_big=0, cape_cnt_surface_msgs=0; - int cape_cnt_no_levels=0, cape_cnt_missing_values=0, cape_cnt_zero_values=0; + int cape_level, cape_count, cape_cnt_too_big, cape_cnt_surface_msgs; + int cape_cnt_no_levels, cape_cnt_missing_values, cape_cnt_zero_values; + int mlcape_count, mlcape_cnt_too_big; + int mlcape_cnt_missing_values, mlcape_cnt_zero_values; float cape_p, cape_h; float cape_qm = bad_data_float; @@ -997,8 +1003,9 @@ void process_pbfile(int i_pb) { bool has_pbl_data; bool do_pbl = false; - bool cal_cape = bufr_obs_name_arr.has(derived_cape, cape_code, false); bool cal_pbl = bufr_obs_name_arr.has(derived_pbl, pbl_code, false); + bool cal_cape = bufr_obs_name_arr.has(derived_cape, cape_code, false); + bool cal_mlcape = bufr_obs_name_arr.has(derived_mlcape, mlcape_code, false); bool is_same_header; unixtime prev_hdr_vld_ut = (unixtime) 0; @@ -1010,6 +1017,10 @@ void process_pbfile(int i_pb) { // Initialize prev_hdr_lat = prev_hdr_lon = prev_hdr_elv = bad_data_double; + cape_level = cape_count = cape_cnt_too_big = 0; + cape_cnt_no_levels = cape_cnt_missing_values = cape_cnt_zero_values = 0; + mlcape_count = mlcape_cnt_too_big = 0; + mlcape_cnt_missing_values = mlcape_cnt_zero_values = 0; if (cal_pbl) { is_same_header = false; @@ -1327,7 +1338,7 @@ void process_pbfile(int i_pb) { } } - if (cal_cape) { + if (cal_cape || cal_mlcape) { cape_level = 0; } @@ -1471,7 +1482,7 @@ void process_pbfile(int i_pb) { obs_arr[4] += c_to_k; } - if (cal_cape) { + if (cal_cape || cal_mlcape) { if(grib_code == pres_grib_code) { is_cape_input = true; if (cape_level < MAX_CAPE_LEVEL) cape_data_pres[cape_level] = obs_arr[4]; @@ -1567,7 +1578,7 @@ void process_pbfile(int i_pb) { } // end for i } - if (cal_cape) { + if (cal_cape || cal_mlcape) { if (cape_member_cnt >= 3) cape_level++; } @@ -1603,7 +1614,7 @@ void process_pbfile(int i_pb) { } } // end for lv - if (cal_cape) { + if (cal_cape || cal_mlcape) { if (1 < cape_level) { bool reverse_levels; float cape_val, cin_val, PLCL,PEQL; @@ -1646,14 +1657,6 @@ void process_pbfile(int i_pb) { } } - //p1d = cape_p; - //t1d = cape_data_temp[cape_level-1]; - //q1d = cape_data_spfh[cape_level-1]; - calcape_(&ivirt,&itype, cape_data_temp, cape_data_spfh, cape_data_pres, - &p1d,&t1d,&q1d, static_dummy_201, - &cape_level, &IMM,&JMM, &cape_level, - &cape_val, &cin_val, &PLCL, &PEQL, static_dummy_200); - if(mlog.verbosity_level() >= 7) { mlog << Debug(7) << method_name << " index,P,T,Q to compute CAPE from " << i_read << "-th message\n" ; @@ -1661,37 +1664,83 @@ void process_pbfile(int i_pb) { mlog << Debug(7) << method_name << " " << idx << ", " << cape_data_pres[idx] << ", " << cape_data_temp[idx] << ", " << cape_data_spfh[idx] << "\n"; } - mlog << Debug(7) << method_name - << " calcape_(" << ivirt << "," << itype << ") cape_val: " - << cape_val << " cape_level: " << cape_level - << ", cin_val: " << cin_val - << ", PLCL: " << PLCL << ", PEQL: " << PEQL - << " lat: " << hdr_lat << ", lon: " << hdr_lon - << " valid_time: " << unix_to_yyyymmdd_hhmmss(hdr_vld_ut) - << " " << hdr_typ << " " << hdr_sid - << "\n\n" ; } - if (cape_val > MAX_CAPE_VALUE) { - cape_cnt_too_big++; - mlog << Debug(5) << method_name - << " Ignored cape_value: " << cape_val << " cape_level: " << cape_level - << ", cin_val: " << cin_val - << ", PLCL: " << PLCL << ", PEQL: " << PEQL << "\n"; + if (cal_cape) { + ivirt = 1; + itype = 1; + calcape_(&ivirt,&itype, cape_data_temp, cape_data_spfh, cape_data_pres, + &p1d,&t1d,&q1d, static_dummy_201, + &cape_level, &IMM,&JMM, &cape_level, + &cape_val, &cin_val, &PLCL, &PEQL, static_dummy_200); + + if(mlog.verbosity_level() >= 7) { + mlog << Debug(7) << method_name + << " calcape_(" << ivirt << "," << itype << ") cape_val: " + << cape_val << " cape_level: " << cape_level + << ", cin_val: " << cin_val + << ", PLCL: " << PLCL << ", PEQL: " << PEQL + << " lat: " << hdr_lat << ", lon: " << hdr_lon + << " valid_time: " << unix_to_yyyymmdd_hhmmss(hdr_vld_ut) + << " " << hdr_typ << " " << hdr_sid + << "\n\n" ; + } + + if (cape_val > MAX_CAPE_VALUE) { + cape_cnt_too_big++; + mlog << Debug(5) << method_name + << " Ignored cape_value: " << cape_val << " cape_level: " << cape_level + << ", cin_val: " << cin_val + << ", PLCL: " << PLCL << ", PEQL: " << PEQL << "\n"; + } + else if (cape_val >= 0) { + obs_arr[1] = cape_code; + obs_arr[2] = cape_p; + obs_arr[3] = cape_h; + obs_arr[4] = cape_val; // observation value + addObservation(obs_arr, (string)hdr_typ, (string)hdr_sid, hdr_vld_ut, + hdr_lat, hdr_lon, hdr_elv, cape_qm, + OBS_BUFFER_SIZE); + cape_count++; + n_derived_obs++; + if (is_eq(cape_val, 0.)) cape_cnt_zero_values++; + } + else cape_cnt_missing_values++; } - else if (cape_val >= 0) { - obs_arr[1] = cape_code; - obs_arr[2] = cape_p; - obs_arr[3] = cape_h; - obs_arr[4] = cape_val; // observation value - addObservation(obs_arr, (string)hdr_typ, (string)hdr_sid, hdr_vld_ut, - hdr_lat, hdr_lon, hdr_elv, cape_qm, - OBS_BUFFER_SIZE); - cape_count++; - n_derived_obs++; - if (is_eq(cape_val, 0.)) cape_cnt_zero_values++; + + if (cal_mlcape) { + float mlcape_val = bad_data_float; + + ivirt = 1; + itype = 2; + // The second last seems to be better than the average of last two or three + p1d = cape_data_pres[cape_level-2]; + t1d = cape_data_temp[cape_level-2]; + q1d = cape_data_spfh[cape_level-2]; + calcape_(&ivirt,&itype, cape_data_temp, cape_data_spfh, cape_data_pres, + &p1d,&t1d,&q1d, static_dummy_201, + &cape_level, &IMM,&JMM, &cape_level, + &mlcape_val, &cin_val, &PLCL, &PEQL, static_dummy_200); + if (mlcape_val > MAX_CAPE_VALUE) { + mlcape_cnt_too_big++; + mlog << Debug(5) << method_name + << " Ignored ML_cape: " << mlcape_val << " cape_level: " + << cape_level << "\n"; + } + else if (mlcape_val >= 0) { + obs_arr[1] = mlcape_code; + obs_arr[2] = cape_p; + obs_arr[3] = cape_h; + obs_arr[4] = mlcape_val; // observation value + addObservation(obs_arr, (string)hdr_typ, (string)hdr_sid, hdr_vld_ut, + hdr_lat, hdr_lon, hdr_elv, cape_qm, + OBS_BUFFER_SIZE); + mlcape_count++; + n_derived_obs++; + if (is_eq(mlcape_val, 0.)) mlcape_cnt_zero_values++; + } + else mlcape_cnt_missing_values++; } - else cape_cnt_missing_values++; } else if (1 < buf_nlev) cape_cnt_no_levels++; else cape_cnt_surface_msgs++; @@ -1948,13 +1997,13 @@ void process_pbfile(int i_pb) { << "Total observations retained or derived\t= " << (n_file_obs + n_derived_obs) << "\n"; - if (cal_cape) { - mlog << Debug(3) << "\nDerived CAPE = " << cape_count - << "\tZero = " << cape_cnt_zero_values + if (cal_cape || cal_mlcape) { + mlog << Debug(3) << "\nDerived CAPE = " << (cape_count + mlcape_count) + << "\tZero = " << (cape_cnt_zero_values + mlcape_cnt_zero_values) << "\n\tnot derived: No cape inputs = " << (cape_cnt_no_levels) << "\tNo vertical levels = " << (cape_cnt_surface_msgs) - << "\n\tfiltered: " << cape_cnt_missing_values << ", " - << cape_cnt_too_big + << "\n\tfiltered: " << (cape_cnt_missing_values + mlcape_cnt_missing_values) << ", " + << (cape_cnt_too_big + mlcape_cnt_too_big) << "\n"; } @@ -3003,7 +3052,11 @@ float compute_pbl(map pqtzuv_map_tq, if (pbl_level <= 0) { mlog << Debug(4) << method_name - << "Skip CALPBL because an empty list after combining TQZ and UV\n"; + << "Skip CALPBL because of an empty list after combining TQZ and UV\n"; + } + else if (pbl_level == 1) { + mlog << Debug(4) << method_name + << "Skip CALPBL because of only one available record after combining TQZ and UV\n"; } else { // Order all observations by pressure from bottom to top @@ -3057,7 +3110,7 @@ float compute_pbl(map pqtzuv_map_tq, pbl_data_hgt[index] = pqtzuv[3]; pbl_data_ugrd[index] = pqtzuv[4]; pbl_data_vgrd[index] = pqtzuv[5]; - mlog << Debug(5) << method_name << "Force to add " + mlog << Debug(6) << method_name << "Force to add " << pres_level << " into " << index << "\n"; index--; } @@ -3077,11 +3130,11 @@ float compute_pbl(map pqtzuv_map_tq, if (hgt_cnt < pbl_level) { hgt_cnt += interpolate_by_pressure(pbl_level, pbl_data_pres, pbl_data_hgt); - mlog << Debug(4) << method_name << "interpolate Z (HGT)\n"; + mlog << Debug(6) << method_name << "interpolate Z (HGT)\n"; } if (spfh_cnt < pbl_level) { spfh_cnt += interpolate_by_pressure(pbl_level, pbl_data_pres, pbl_data_spfh); - mlog << Debug(4) << method_name << "interpolate Q (SPFH)\n"; + mlog << Debug(6) << method_name << "interpolate Q (SPFH)\n"; } if ((spfh_cnt>0) && (pbl_level>0)) { diff --git a/met/src/tools/other/pb2nc/pb2nc_conf_info.cc b/met/src/tools/other/pb2nc/pb2nc_conf_info.cc index d940939072..3c3a53c74a 100644 --- a/met/src/tools/other/pb2nc/pb2nc_conf_info.cc +++ b/met/src/tools/other/pb2nc/pb2nc_conf_info.cc @@ -25,6 +25,18 @@ using namespace std; #include "grib_strings.h" #include "vx_log.h" + +//////////////////////////////////////////////////////////////////////// + +map parse_conf_obs_bufr_map(Dictionary *dict) { + + const char *key_name = (0 != dict->lookup_array(conf_key_obs_prepbufr_map, false, false)) + ? conf_key_obs_prepbufr_map : conf_key_old_prepbufr_map; + map m = parse_conf_key_value_map(dict, (const char *)key_name); + parse_add_conf_key_value_map(dict, conf_key_obs_bufr_map, &m); + return m; +} + //////////////////////////////////////////////////////////////////////// // // Code for class PB2NCConfInfo @@ -85,14 +97,42 @@ void PB2NCConfInfo::clear() { void PB2NCConfInfo::read_config(const char *default_file_name, const char *user_file_name) { + int warning_code = 0; + bool use_bad_one = true; + string bad_file_names; + // Read the config file constants conf.read(replace_path(config_const_filename).c_str()); // Read the default config file conf.read(default_file_name); + if (0 != conf.lookup_array(conf_key_old_prepbufr_map, false, false)) { + warning_code = 1; + bad_file_names = default_file_name; + if (0 != conf.lookup_array(conf_key_obs_prepbufr_map, false, false)) use_bad_one = false; + } // Read the user-specified config file conf.read(user_file_name); + if (0 != conf.lookup_array(conf_key_old_prepbufr_map, false, false)) { + warning_code = 2; + if (0 < bad_file_names.length()) bad_file_names += " and "; + bad_file_names += user_file_name; + if (0 != conf.lookup_array(conf_key_obs_prepbufr_map, false, false)) use_bad_one = false; + } + + if (0 < warning_code) { + string ignored_message = " "; + + if (!use_bad_one) { + ignored_message = " (Ignored \""; + ignored_message += conf_key_old_prepbufr_map; + ignored_message += "\" key)"; + } + mlog << Warning << "\nPlease rename the configuration key \"" + << conf_key_old_prepbufr_map << "\" to \"" << conf_key_obs_prepbufr_map + << "\" at " << bad_file_names << ignored_message << "\n\n"; + } return; } diff --git a/test/config/PB2NCConfig_airnow b/test/config/PB2NCConfig_airnow index f9cba2c0bb..5fd9a2b01d 100644 --- a/test/config/PB2NCConfig_airnow +++ b/test/config/PB2NCConfig_airnow @@ -94,7 +94,7 @@ obs_bufr_map = []; // This map is for PREPBUFR. It will be added into obs_bufr_map. // Please do not override this map. // -obs_prefbufr_map = [ +obs_prepbufr_map = [ { key = "POB"; val = "PRES"; }, { key = "QOB"; val = "SPFH"; }, { key = "TOB"; val = "TMP"; }, diff --git a/test/config/PB2NCConfig_pbl b/test/config/PB2NCConfig_pbl index eeedd7a3e4..ac48b5850f 100644 --- a/test/config/PB2NCConfig_pbl +++ b/test/config/PB2NCConfig_pbl @@ -92,7 +92,7 @@ level_category = [0, 1, 4, 5, 6]; // Use obs_bufr_map to rename variables in the output. // If empty, process all available variables. // -obs_bufr_var = ["D_CAPE", "D_PBL"]; +obs_bufr_var = ["D_CAPE", "D_PBL", "D_MLCAPE"]; //////////////////////////////////////////////////////////////////////////////// // @@ -106,7 +106,7 @@ obs_bufr_map = []; // abbreviations in the output. This default map is appended to obs_bufr_map. // This should not typically be overridden. // -obs_prefbufr_map = [ +obs_prepbufr_map = [ { key = "POB"; val = "PRES"; }, { key = "QOB"; val = "SPFH"; }, { key = "TOB"; val = "TMP"; }, @@ -120,6 +120,7 @@ obs_prefbufr_map = [ { key = "D_PBL"; val = "HPBL"; }, { key = "D_PRMSL"; val = "PRMSL"; }, { key = "D_CAPE"; val = "CAPE"; }, + { key = "D_MLCAPE"; val = "MLCAPE"; }, { key = "TDO"; val = "DPT"; }, { key = "PMO"; val = "PRMSL"; }, { key = "TOCC"; val = "TCDC"; },