PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
OrographicPrecipitationSerial.cc
Go to the documentation of this file.
1 // Copyright (C) 2018, 2019, 2020, 2021, 2023 Andy Aschwanden and Constantine Khroulev
2 //
3 // This file is part of PISM.
4 //
5 // PISM is free software; you can redistribute it and/or modify it under the
6 // terms of the GNU General Public License as published by the Free Software
7 // Foundation; either version 3 of the License, or (at your option) any later
8 // version.
9 //
10 // PISM is distributed in the hope that it will be useful, but WITHOUT ANY
11 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 // details.
14 //
15 // You should have received a copy of the GNU General Public License
16 // along with PISM; if not, write to the Free Software
17 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18 
19 #include "pism/coupler/atmosphere/OrographicPrecipitationSerial.hh"
20 
21 #include <complex> // std::complex<double>, std::sqrt()
22 #include <gsl/gsl_math.h> // M_PI
23 #include <cassert> // assert()
24 #include <cmath> // std::exp()
25 
26 #include "pism/util/ConfigInterface.hh"
27 #include "pism/util/error_handling.hh"
28 #include "pism/util/fftw_utilities.hh"
29 
30 namespace pism {
31 namespace atmosphere {
32 
33 /*!
34  * @param[in] config configuration database
35  * @param[in] log logger
36  * @param[in] Mx grid size in the X direction
37  * @param[in] My grid size in the Y direction
38  * @param[in] dx grid spacing in the X direction
39  * @param[in] dy grid spacing in the Y direction
40  * @param[in] Nx extended grid size in the X direction
41  * @param[in] Ny extended grid size in the Y direction
42  */
44  int Mx, int My,
45  double dx, double dy,
46  int Nx, int Ny)
47  : m_Mx(Mx), m_My(My), m_Nx(Nx), m_Ny(Ny) {
48 
49  m_eps = 1.0e-18;
50 
51  // derive more parameters
52  {
53  m_i0_offset = (Nx - Mx) / 2;
54  m_j0_offset = (Ny - My) / 2;
55 
56  m_kx = fftfreq(m_Nx, dx / (2.0 * M_PI));
57  m_ky = fftfreq(m_Ny, dy / (2.0 * M_PI));
58  }
59 
60  {
61  m_background_precip_pre = config.get_number("atmosphere.orographic_precipitation.background_precip_pre", "mm/s");
62  m_background_precip_post = config.get_number("atmosphere.orographic_precipitation.background_precip_post", "mm/s");
63 
64  m_precip_scale_factor = config.get_number("atmosphere.orographic_precipitation.scale_factor");
65  m_tau_c = config.get_number("atmosphere.orographic_precipitation.conversion_time");
66  m_tau_f = config.get_number("atmosphere.orographic_precipitation.fallout_time");
67  m_Hw = config.get_number("atmosphere.orographic_precipitation.water_vapor_scale_height");
68  m_Nm = config.get_number("atmosphere.orographic_precipitation.moist_stability_frequency");
69  m_wind_speed = config.get_number("atmosphere.orographic_precipitation.wind_speed");
70  m_wind_direction = config.get_number("atmosphere.orographic_precipitation.wind_direction");
71  m_gamma = config.get_number("atmosphere.orographic_precipitation.lapse_rate");
72  m_Theta_m = config.get_number("atmosphere.orographic_precipitation.moist_adiabatic_lapse_rate");
73  m_rho_Sref = config.get_number("atmosphere.orographic_precipitation.reference_density");
74  m_latitude = config.get_number("atmosphere.orographic_precipitation.coriolis_latitude");
75  m_truncate = config.get_flag("atmosphere.orographic_precipitation.truncate");
76 
77 
78  // derived constants
79  m_f = 2.0 * 7.2921e-5 * sin(m_latitude * M_PI / 180.0);
80 
81  m_u = -sin(m_wind_direction * 2.0 * M_PI / 360.0) * m_wind_speed;
82  m_v = -cos(m_wind_direction * 2.0 * M_PI / 360.0) * m_wind_speed;
83 
85  }
86 
87  // memory allocation
88  {
89  PetscErrorCode ierr = 0;
90 
91  // precipitation
92  ierr = VecCreateSeq(PETSC_COMM_SELF, m_Mx * m_My, m_precipitation.rawptr());
93  PISM_CHK(ierr, "VecCreateSeq");
94 
95  // FFTW arrays
96  m_fftw_input = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * m_Nx * m_Ny);
97  m_fftw_output = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * m_Nx * m_Ny);
98  m_G_hat = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * m_Nx * m_Ny);
99 
100  // FFTW plans
101  m_dft_forward = fftw_plan_dft_2d(m_Nx, m_Ny, m_fftw_input, m_fftw_output,
102  FFTW_FORWARD, FFTW_ESTIMATE);
103  m_dft_inverse = fftw_plan_dft_2d(m_Nx, m_Ny, m_fftw_input, m_fftw_output,
104  FFTW_BACKWARD, FFTW_ESTIMATE);
105 
106  // Note: FFTW is weird. If a malloc() call fails it will just call
107  // abort() on you without giving you a chance to recover or tell the
108  // user what happened. This is why we don't check return values of
109  // fftw_malloc() and fftw_plan_dft_2d() calls here...
110  //
111  // (Constantine Khroulev, February 1, 2015)
112  }
113 
114  // initialize the Gaussian filter
115  {
116  FFTWArray G_hat(m_G_hat, m_Nx, m_Ny);
117  double sigma = config.get_number("atmosphere.orographic_precipitation.smoothing_standard_deviation");
118 
119  if (sigma > 0.0) {
120  FFTWArray
121  fftw_output(m_fftw_output, m_Nx, m_Ny),
122  fftw_input(m_fftw_input, m_Nx, m_Ny);
123 
124  int
125  Nx2 = Nx / 2,
126  Ny2 = Ny / 2;
127 
128  double sum = 0.0;
129  for (int i = 0; i < m_Nx; i++) {
130  for (int j = 0; j < m_Ny; j++) {
131  int
132  p = i <= Nx2 ? i : m_Nx - i,
133  q = j <= Ny2 ? j : m_Ny - j;
134  double
135  x = p * dx,
136  y = q * dy;
137 
138  double G = std::exp(-0.5 * (x * x + y * y) / (sigma * sigma));
139  sum += G;
140 
141  fftw_input(i, j) = G;
142  }
143  }
144 
145  // normalize:
146  assert(sum > 0.0);
147  for (int i = 0; i < m_Nx; i++) {
148  for (int j = 0; j < m_Ny; j++) {
149  fftw_input(i, j) /= sum;
150  }
151  }
152 
153  // compute FFT of the Gaussian
154  fftw_execute(m_dft_forward);
155 
156  // copy to m_G_hat
157  for (int i = 0; i < m_Nx; i++) {
158  for (int j = 0; j < m_Ny; j++) {
159  G_hat(i, j) = fftw_output(i, j);
160  }
161  }
162 
163  } else {
164  // fill m_G_hat with ones to disable smoothing
165  for (int i = 0; i < m_Nx; i++) {
166  for (int j = 0; j < m_Ny; j++) {
167  G_hat(i, j) = 1.0;
168  }
169  }
170  }
171  }
172 }
173 
175  fftw_destroy_plan(m_dft_forward);
176  fftw_destroy_plan(m_dft_inverse);
177  fftw_free(m_fftw_input);
178  fftw_free(m_fftw_output);
179  fftw_free(m_G_hat);
180 }
181 
182 /*!
183  * Return precipitation (FIXME: units?)
184  */
186  return m_precipitation;
187 }
188 
189 /*!
190  * Update precipitation.
191  *
192  * @param[in] surface_elevation surface on the physical (Mx*My) grid
193  */
195  // solves:
196  // Phat(k,l) = (Cw * i * sigma * Hhat(k,l)) /
197  // (1 - i * m * Hw) * (1 + i * sigma * tauc) * (1 + i * sigma * tauf);
198  // see equation (49) in
199  // R. B. Smith and I. Barstad, 2004:
200  // A Linear Theory of Orographic Precipitation. J. Atmos. Sci. 61, 1377-1391.
201 
202  std::complex<double> I(0.0, 1.0);
203 
204  // Compute fft2(surface_elevation)
205  {
207  set_real_part(surface_elevation,
208  1.0,
209  m_Mx, m_My,
210  m_Nx, m_Ny,
212  m_fftw_input);
213  fftw_execute(m_dft_forward);
214  }
215 
216  {
217  FFTWArray
218  fftw_output(m_fftw_output, m_Nx, m_Ny),
219  fftw_input(m_fftw_input, m_Nx, m_Ny),
220  G_hat(m_G_hat, m_Nx, m_Ny);
221 
222  for (int i = 0; i < m_Nx; i++) {
223  const double kx = m_kx[i];
224  for (int j = 0; j < m_Ny; j++) {
225  const double ky = m_ky[j];
226 
227  // FFT(h) * FFT(Gaussian), i.e. FFT(smoothed ice surface elevation)
228  const auto &h_hat = fftw_output(i, j) * G_hat(i, j);
229 
230  double sigma = m_u * kx + m_v * ky;
231 
232  // See equation (6) in [@ref SmithBarstadBonneau2005]
233  std::complex<double> m;
234  {
235  double denominator = sigma * sigma - m_f * m_f;
236 
237  // avoid dividing by zero:
238  if (fabs(denominator) < m_eps) {
239  denominator = denominator >= 0 ? m_eps : -m_eps;
240  }
241 
242  double m_squared = (m_Nm * m_Nm - sigma * sigma) * (kx * kx + ky * ky) / denominator;
243 
244  // Note: this is a *complex* square root.
245  m = std::sqrt(std::complex<double>(m_squared));
246 
247  if (m_squared >= 0.0 and sigma != 0.0) {
248  m *= sigma > 0.0 ? 1.0 : -1.0;
249  }
250  }
251 
252  // avoid dividing by zero:
253  double delta = 0.0;
254  if (std::abs(1.0 - I * m * m_Hw) < m_eps) {
255  delta = m_eps;
256  }
257 
258  // See equation (49) in [@ref SmithBarstad2004] or equation (3) in [@ref
259  // SmithBarstadBonneau2005].
260  auto P_hat = h_hat * (m_Cw * I * sigma /
261  ((1.0 - I * m * m_Hw + delta) *
262  (1.0 + I * sigma * m_tau_c) *
263  (1.0 + I * sigma * m_tau_f)));
264  // Note: sigma, m_tau_c, and m_tau_f are purely real, so the second and the third
265  // factors in the denominator are never zero.
266  //
267  // The first factor (1 - i m H_w) *could* be zero. Here we check if it is and
268  // "regularize" if necessary.
269 
270  fftw_input(i, j) = P_hat;
271  }
272  }
273  }
274 
275  fftw_execute(m_dft_inverse);
276 
277  // get m_fftw_output and put it into m_precipitation
279  1.0 / (m_Nx * m_Ny),
280  m_Mx, m_My,
281  m_Nx, m_Ny,
284 
286  for (int i = 0; i < m_Mx; i++) {
287  for (int j = 0; j < m_My; j++) {
288  p(i, j) += m_background_precip_pre;
289  if (m_truncate) {
290  p(i, j) = std::max(p(i, j), 0.0);
291  }
292  p(i, j) *= m_precip_scale_factor;
293  p(i, j) += m_background_precip_post;
294  }
295  }
296 }
297 
298 } // end of namespace atmosphere
299 } // end of namespace pism
double get_number(const std::string &name, UseFlag flag=REMEMBER_THIS_USE) const
bool get_flag(const std::string &name, UseFlag flag=REMEMBER_THIS_USE) const
A class for storing and accessing PISM configuration flags and parameters.
T * rawptr()
Definition: Wrapper.hh:39
OrographicPrecipitationSerial(const Config &config, int Mx, int My, double dx, double dy, int Nx, int Ny)
Wrapper around VecGetArray2d and VecRestoreArray2d.
Definition: Vec.hh:55
#define PISM_CHK(errcode, name)
const double G
Definition: exactTestK.c:40
double max(const array::Scalar &input)
Finds maximum over all the values in an array::Scalar object. Ignores ghosts.
Definition: Scalar.cc:165
double sum(const array::Scalar &input)
Sums up all the values in an array::Scalar object. Ignores ghosts.
Definition: Scalar.cc:150
void clear_fftw_array(fftw_complex *input, int Nx, int Ny)
Fill input with zeros.
void get_real_part(fftw_complex *input, double normalization, int Mx, int My, int Nx, int Ny, int i0, int j0, petsc::Vec &output)
Get the real part of input and put it in output.
void set_real_part(petsc::Vec &input, double normalization, int Mx, int My, int Nx, int Ny, int i0, int j0, fftw_complex *output)
std::vector< double > fftfreq(int M, double normalization)
const int I[]
Definition: ssafd_code.cc:24