Loading [MathJax]/extensions/tex2jax.js
PISM, A Parallel Ice Sheet Model 2.2.2-d6b3a29ca committed by Constantine Khrulev on 2025-03-28
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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
30namespace pism {
31namespace 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) {
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 */
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,
213 fftw_execute(m_dft_forward);
214 }
215
216 {
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
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)