UFO
MWCLWCheck.cc
Go to the documentation of this file.
1 /*
2  * (C) Copyright 2018-2019 UCAR
3  *
4  * This software is licensed under the terms of the Apache Licence Version 2.0
5  * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
6  */
7 
9 
10 #include <math.h>
11 
12 #include <algorithm>
13 #include <iostream>
14 #include <vector>
15 
16 #include "eckit/config/Configuration.h"
17 
18 #include "ioda/ObsDataVector.h"
19 #include "ioda/ObsSpace.h"
20 #include "oops/util/Logger.h"
21 #include "oops/util/missingValues.h"
22 #include "ufo/filters/QCflags.h"
23 
24 namespace ufo {
25 
26 // -----------------------------------------------------------------------------
27 
28 MWCLWCheck::MWCLWCheck(ioda::ObsSpace & obsdb, const eckit::Configuration & config,
29  std::shared_ptr<ioda::ObsDataVector<int> > flags,
30  std::shared_ptr<ioda::ObsDataVector<float> > obserr)
31  : FilterBase(obsdb, config, flags, obserr), invars_(config_, "clw variables") {
32  oops::Log::debug() << "MWCLWCheck: config = " << config_ << std::endl;
33  const Variable var0(invars_[0] + "@HofX");
34  const Variable var1(invars_[1] + "@HofX");
35  allvars_ += var0;
36  allvars_ += var1;
37 }
38 
39 // -----------------------------------------------------------------------------
40 
42 
43 // -----------------------------------------------------------------------------
44 
45 void MWCLWCheck::applyFilter(const std::vector<bool> & apply,
46  const Variables & filtervars,
47  std::vector<std::vector<bool>> & flagged) const {
48  oops::Log::trace() << "MWCLWCheck postFilter" << std::endl;
49 
50  const oops::Variables observed = obsdb_.obsvariables();
51  const float missing = util::missingValue(missing);
52  float amsua_clw(float, float, float);
53 
54 // Get config
55  std::vector<float> clw_thresholds = config_.getFloatVector("clw_thresholds");
56 // clw_option controls how the clw is calculated:
57 // 1) Use observed BTs
58 // 2) Use calculated BTs
59 // 3) Symmetric calculation
60  const int clw_option = config_.getInt("clw_option");
61  oops::Log::debug() << "MWCLWCheck: config = " << config_ << std::endl;
62 
63 // Number of channels to be tested and number of thresholds needs to be the same
64  ASSERT(clw_thresholds.size() == filtervars.nvars());
65 // Thresholds need to have non-missing values
66  for (size_t jv = 0; jv < filtervars.nvars(); ++jv) {
67  ASSERT(clw_thresholds[jv] != missing);
68  }
69 // Check we have the correct number of channels to do the CLW calculation
70  ASSERT(invars_.size() == 2);
71 // Check clw_option is in range
72  ASSERT(clw_option >= 1 && clw_option <=3);
73 
74  ioda::ObsDataVector<float> obs(obsdb_, filtervars.toOopsVariables(), "ObsValue");
75  ioda::ObsDataVector<float> obs_for_calc(obsdb_, invars_, "ObsValue");
76  ioda::ObsDataVector<float> sza(obsdb_, "sensor_zenith_angle", "MetaData");
77  ioda::ObsDataVector<float> clw(obsdb_, "cloud_liquid_water", "Diagnostic", false);
78  ioda::ObsDataVector<float> clw_guess_out(obsdb_, "clws_guess", "Diagnostic", false);
79  ioda::ObsDataVector<float> clw_obs_out(obsdb_, "clw_obs", "Diagnostic", false);
80 
81 // H(x)
82  const Variable var0(invars_[0] + "@HofX");
83  const Variable var1(invars_[1] + "@HofX");
84  std::vector<float> hofx0;
85  data_.get(var0, hofx0);
86  std::vector<float> hofx1;
87  data_.get(var1, hofx1);
88 
89  auto clw_obs = [&](size_t jobs) { return amsua_clw(
90  obs_for_calc[0][jobs], obs_for_calc[1][jobs], sza[0][jobs]); };
91 
92  auto clw_guess = [&](size_t jobs) { return amsua_clw(
93  hofx0[jobs], hofx1[jobs], sza[0][jobs]); };
94 
95 // Loop over obs locations calculating CLW from observations
96  for (size_t jobs = 0; jobs < obs.nlocs(); ++jobs) {
97  switch (clw_option) {
98  case 1 :
99  clw[0][jobs] = clw_obs(jobs);
100  break;
101  case 2 :
102  clw[0][jobs] = clw_guess(jobs);
103  break;
104  case 3 :
105  clw_obs_out[0][jobs] = clw_obs(jobs);
106  clw_guess_out[0][jobs] = clw_guess(jobs);
107  if (clw_obs_out[0][jobs] == missing || clw_guess_out[0][jobs] == missing) {
108  clw[0][jobs] = missing;
109  } else {
110  clw[0][jobs] = (clw_obs_out[0][jobs] + clw_guess_out[0][jobs])/2.0;
111  }
112  break;
113  default:
114  oops::Log::error() << "Invalid value for clw_option:" << clw_option << std::endl;
115  ABORT("Invalid value for clw_option");
116  }
117 
118 // Apply CLW threshold to observations
119  for (size_t jv = 0; jv < filtervars.nvars(); ++jv) {
120  if (apply[jobs]) {
121  if (clw[0][jobs] == missing ||
122  clw[0][jobs] > clw_thresholds[jv]) flagged[jv][jobs] = true;
123  }
124  }
125  }
126  clw.save("Derived");
127  clw_obs_out.save("Derived");
128  clw_guess_out.save("Derived");
129 }
130 
131 // -----------------------------------------------------------------------------
132 
133 float amsua_clw(float tobs1, float tobs2, float sza) {
134  const float d1 = 0.754;
135  const float d2 = -2.265;
136  const float missing = util::missingValue(missing);
137 
138  float clw;
139 
140  if (tobs1 != missing && tobs2 != missing && sza != missing &&
141  tobs1 <= 284.0 && tobs2 <= 284.0 && tobs1 > 0.0 && tobs2 > 0.0) {
142  float cossza = cos(M_PI * sza/180.0);
143  float d0 = 8.240 - (2.622 - 1.846*cossza)*cossza;
144  clw = cossza*(d0 + d1*log(285.0-tobs1)) + d2*log(285.0-tobs2);
145  clw = std::max(0.0f, clw);
146  } else {
147  clw = missing;
148  }
149 
150  return clw;
151 }
152 
153 // -----------------------------------------------------------------------------
154 
155 void MWCLWCheck::print(std::ostream & os) const {
156  os << "MWCLWCheck: config = " << config_ << std::endl;
157 }
158 
159 // -----------------------------------------------------------------------------
160 
161 } // namespace ufo
ufo::Variables::nvars
size_t nvars() const
Definition: Variables.cc:104
ufo::Variables
Definition: src/ufo/filters/Variables.h:24
ufo::MWCLWCheck::MWCLWCheck
MWCLWCheck(ioda::ObsSpace &, const eckit::Configuration &, std::shared_ptr< ioda::ObsDataVector< int > >, std::shared_ptr< ioda::ObsDataVector< float > >)
Definition: MWCLWCheck.cc:28
ufo_radiancerttov_utils_mod::debug
logical, public debug
Definition: ufo_radiancerttov_utils_mod.F90:100
ufo::amsua_clw
float amsua_clw(float tobs1, float tobs2, float sza)
Definition: MWCLWCheck.cc:133
ufo::FilterBase::obsdb_
ioda::ObsSpace & obsdb_
Definition: FilterBase.h:59
ufo::FilterBase
FilterBase: Base class for UFO QC filters.
Definition: FilterBase.h:42
ufo
Definition: RunCRTM.h:27
ufo::FilterBase::data_
ObsFilterData data_
Definition: FilterBase.h:65
ufo::MWCLWCheck::invars_
oops::Variables invars_
Definition: MWCLWCheck.h:52
ufo::MWCLWCheck::print
void print(std::ostream &) const override
Definition: MWCLWCheck.cc:155
ufo::FilterBase::allvars_
ufo::Variables allvars_
Definition: FilterBase.h:63
ufo::MWCLWCheck::applyFilter
void applyFilter(const std::vector< bool > &, const Variables &, std::vector< std::vector< bool >> &) const override
Definition: MWCLWCheck.cc:45
ufo::QCflags::missing
constexpr int missing
Definition: QCflags.h:15
MWCLWCheck.h
ufo::Variables::toOopsVariables
oops::Variables toOopsVariables() const
Definition: Variables.cc:144
QCflags.h
ufo::QCflags::clw
constexpr int clw
Definition: QCflags.h:23
ioda::ObsDataVector< int >
ufo::FilterBase::config_
const eckit::LocalConfiguration config_
Definition: FilterBase.h:60
ufo::MWCLWCheck::~MWCLWCheck
~MWCLWCheck()
Definition: MWCLWCheck.cc:41
ufo::ObsFilterData::get
void get(const Variable &, std::vector< float > &) const
Gets requested data from ObsFilterData.
Definition: ObsFilterData.cc:130
ufo::Variable
Definition: Variable.h:23