tOrographicPrecipitationSerial.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tOrographicPrecipitationSerial.cc (8596B)
---
1 // Copyright (C) 2018, 2019, 2020 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 "OrographicPrecipitationSerial.hh"
20
21 #include <complex> // std::complex<double>, std::sqrt()
22 #include <fftw3.h>
23 #include <gsl/gsl_math.h> // M_PI
24
25 #include "pism/util/ConfigInterface.hh"
26 #include "pism/util/error_handling.hh"
27 #include "pism/util/petscwrappers/Vec.hh"
28 #include "pism/util/pism_utilities.hh"
29 #include "pism/util/fftw_utilities.hh"
30
31 namespace pism {
32 namespace atmosphere {
33
34 /*!
35 * @param[in] config configuration database
36 * @param[in] log logger
37 * @param[in] Mx grid size in the X direction
38 * @param[in] My grid size in the Y direction
39 * @param[in] dx grid spacing in the X direction
40 * @param[in] dy grid spacing in the Y direction
41 * @param[in] Nx extended grid size in the X direction
42 * @param[in] Ny extended grid size in the Y direction
43 */
44 OrographicPrecipitationSerial::OrographicPrecipitationSerial(const Config &config,
45 int Mx, int My,
46 double dx, double dy,
47 int Nx, int Ny)
48 : m_Mx(Mx), m_My(My), m_Nx(Nx), m_Ny(Ny) {
49
50 m_eps = 1.0e-18;
51
52 // derive more parameters
53 {
54 m_i0_offset = (Nx - Mx) / 2;
55 m_j0_offset = (Ny - My) / 2;
56
57 m_kx = fftfreq(m_Nx, dx / (2.0 * M_PI));
58 m_ky = fftfreq(m_Ny, dy / (2.0 * M_PI));
59 }
60
61 {
62 m_background_precip_pre = config.get_number("atmosphere.orographic_precipitation.background_precip_pre", "mm/s");
63 m_background_precip_post = config.get_number("atmosphere.orographic_precipitation.background_precip_post", "mm/s");
64
65 m_precip_scale_factor = config.get_number("atmosphere.orographic_precipitation.scale_factor");
66 m_tau_c = config.get_number("atmosphere.orographic_precipitation.conversion_time");
67 m_tau_f = config.get_number("atmosphere.orographic_precipitation.fallout_time");
68 m_Hw = config.get_number("atmosphere.orographic_precipitation.water_vapor_scale_height");
69 m_Nm = config.get_number("atmosphere.orographic_precipitation.moist_stability_frequency");
70 m_wind_speed = config.get_number("atmosphere.orographic_precipitation.wind_speed");
71 m_wind_direction = config.get_number("atmosphere.orographic_precipitation.wind_direction");
72 m_gamma = config.get_number("atmosphere.orographic_precipitation.lapse_rate");
73 m_Theta_m = config.get_number("atmosphere.orographic_precipitation.moist_adiabatic_lapse_rate");
74 m_rho_Sref = config.get_number("atmosphere.orographic_precipitation.reference_density");
75 m_latitude = config.get_number("atmosphere.orographic_precipitation.coriolis_latitude");
76 m_truncate = config.get_flag("atmosphere.orographic_precipitation.truncate");
77
78
79 // derived constants
80 m_f = 2.0 * 7.2921e-5 * sin(m_latitude * M_PI / 180.0);
81
82 m_u = -sin(m_wind_direction * 2.0 * M_PI / 360.0) * m_wind_speed;
83 m_v = -cos(m_wind_direction * 2.0 * M_PI / 360.0) * m_wind_speed;
84
85 m_Cw = m_rho_Sref * m_Theta_m / m_gamma;
86 }
87
88 // memory allocation
89 {
90 PetscErrorCode ierr = 0;
91
92 // precipitation
93 ierr = VecCreateSeq(PETSC_COMM_SELF, m_Mx * m_My, m_precipitation.rawptr());
94 PISM_CHK(ierr, "VecCreateSeq");
95
96 // FFTW arrays
97 m_fftw_input = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * m_Nx * m_Ny);
98 m_fftw_output = (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
115 OrographicPrecipitationSerial::~OrographicPrecipitationSerial() {
116 fftw_destroy_plan(m_dft_forward);
117 fftw_destroy_plan(m_dft_inverse);
118 fftw_free(m_fftw_input);
119 fftw_free(m_fftw_output);
120 }
121
122 /*!
123 * Return precipitation (FIXME: units?)
124 */
125 Vec OrographicPrecipitationSerial::precipitation() const {
126 return m_precipitation;
127 }
128
129 /*!
130 * Update precipitation.
131 *
132 * @param[in] surface_elevation surface on the physical (Mx*My) grid
133 */
134 void OrographicPrecipitationSerial::update(Vec surface_elevation) {
135 // solves:
136 // Phat(k,l) = (Cw * i * sigma * Hhat(k,l)) /
137 // (1 - i * m * Hw) * (1 + i * sigma * tauc) * (1 + i * sigma * tauc);
138 // see equation (49) in
139 // R. B. Smith and I. Barstad, 2004:
140 // A Linear Theory of Orographic Precipitation. J. Atmos. Sci. 61, 1377-1391.
141
142 std::complex<double> I(0.0, 1.0);
143
144 // Compute fft2(surface_elevation)
145 {
146 clear_fftw_array(m_fftw_input, m_Nx, m_Ny);
147 set_real_part(surface_elevation,
148 1.0,
149 m_Mx, m_My,
150 m_Nx, m_Ny,
151 m_i0_offset, m_j0_offset,
152 m_fftw_input);
153 fftw_execute(m_dft_forward);
154 }
155
156 {
157 FFTWArray
158 fftw_output(m_fftw_output, m_Nx, m_Ny),
159 fftw_input(m_fftw_input, m_Nx, m_Ny);
160
161 for (int i = 0; i < m_Nx; i++) {
162 const double kx = m_kx[i];
163 for (int j = 0; j < m_Ny; j++) {
164 const double ky = m_ky[j];
165
166 const auto &h_hat = fftw_output(i, j);
167
168 double sigma = m_u * kx + m_v * ky;
169
170 // See equation (6) in [@ref SmithBarstadBonneau2005]
171 std::complex<double> m;
172 {
173 double denominator = sigma * sigma - m_f * m_f;
174
175 // avoid dividing by zero:
176 if (fabs(denominator) < m_eps) {
177 denominator = denominator >= 0 ? m_eps : -m_eps;
178 }
179
180 double m_squared = (m_Nm * m_Nm - sigma * sigma) * (kx * kx + ky * ky) / denominator;
181
182 // Note: this is a *complex* square root.
183 m = std::sqrt(std::complex<double>(m_squared));
184
185 if (m_squared >= 0.0 and sigma != 0.0) {
186 m *= sigma > 0.0 ? 1.0 : -1.0;
187 }
188 }
189
190 // avoid dividing by zero:
191 double delta = 0.0;
192 if (std::abs(1.0 - I * m * m_Hw) < m_eps) {
193 delta = m_eps;
194 }
195
196 // See equation (49) in [@ref SmithBarstad2004] or equation (3) in [@ref
197 // SmithBarstadBonneau2005].
198 auto P_hat = h_hat * (m_Cw * I * sigma /
199 ((1.0 - I * m * m_Hw + delta) *
200 (1.0 + I * sigma * m_tau_c) *
201 (1.0 + I * sigma * m_tau_f)));
202 // Note: sigma, m_tau_c, and m_tau_f are purely real, so the second and the third
203 // factors in the denominator are never zero.
204 //
205 // The first factor (1 - i m H_w) *could* be zero. Here we check if it is and
206 // "regularize" if necessary.
207
208 fftw_input(i, j) = P_hat;
209 }
210 }
211 }
212
213 fftw_execute(m_dft_inverse);
214
215 // get m_fftw_output and put it into m_p
216 get_real_part(m_fftw_output,
217 1.0 / (m_Nx * m_Ny),
218 m_Mx, m_My,
219 m_Nx, m_Ny,
220 m_i0_offset, m_j0_offset,
221 m_precipitation);
222
223 petsc::VecArray2D p(m_precipitation, m_Mx, m_My);
224 for (int i = 0; i < m_Mx; i++) {
225 for (int j = 0; j < m_My; j++) {
226 p(i, j) += m_background_precip_pre;
227 if (m_truncate) {
228 p(i, j) = std::max(p(i, j), 0.0);
229 }
230 p(i, j) *= m_precip_scale_factor;
231 p(i, j) += m_background_precip_post;
232 }
233 }
234 }
235
236 } // end of namespace atmosphere
237 } // end of namespace pism