tSIAFD_Regional.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
---
tSIAFD_Regional.cc (3157B)
---
1 /* Copyright (C) 2015, 2016, 2017, 2019 PISM Authors
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
20 #include "SIAFD_Regional.hh"
21 #include "pism/stressbalance/StressBalance.hh"
22 #include "pism/geometry/Geometry.hh"
23
24 namespace pism {
25
26 namespace stressbalance {
27
28 SIAFD_Regional::SIAFD_Regional(IceGrid::ConstPtr grid)
29 : SIAFD(grid),
30 m_h_x_no_model(grid, "h_x_no_model", WITH_GHOSTS),
31 m_h_y_no_model(grid, "h_y_no_model", WITH_GHOSTS) {
32 // empty
33 }
34
35 SIAFD_Regional::~SIAFD_Regional() {
36 // empty
37 }
38
39 void SIAFD_Regional::init() {
40
41 SIAFD::init();
42
43 m_log->message(2, " using the regional version of the SIA solver...\n");
44 }
45
46 void SIAFD_Regional::compute_surface_gradient(const Inputs &inputs,
47 IceModelVec2Stag &h_x, IceModelVec2Stag &h_y) {
48
49 SIAFD::compute_surface_gradient(inputs, h_x, h_y);
50
51 // this call updates ghosts of h_x_no_model and h_y_no_model
52 surface_gradient_haseloff(*inputs.no_model_surface_elevation,
53 inputs.geometry->cell_type,
54 m_h_x_no_model, m_h_y_no_model);
55
56 const IceModelVec2Int &no_model = *inputs.no_model_mask;
57
58 const int Mx = m_grid->Mx(), My = m_grid->My();
59
60 IceModelVec::AccessList list{&h_x, &h_y, &no_model, &m_h_x_no_model, &m_h_y_no_model};
61
62 for (PointsWithGhosts p(*m_grid); p; p.next()) {
63 const int i = p.i(), j = p.j();
64
65 auto M = no_model.int_box(i, j);
66
67 // x-component, i-offset
68 if (M.ij > 0.5 or M.e > 0.5) {
69
70 if (i < 0 or i + 1 > Mx - 1) {
71 h_x(i, j, 0) = 0.0;
72 } else {
73 h_x(i, j, 0) = m_h_x_no_model(i, j, 0);
74 }
75 }
76
77 // x-component, j-offset
78 if (M.nw > 0.5 or M.ne > 0.5 or M.w > 0.5 or M.e > 0.5) {
79
80 if (i - 1 < 0 or j + 1 > My - 1 or i + 1 > Mx - 1) {
81 h_x(i, j, 1) = 0.0;
82 } else {
83 h_x(i, j, 1) = m_h_x_no_model(i, j, 1);
84 }
85
86 }
87
88 // y-component, i-offset
89 if (M.n > 0.5 or M.ne > 0.5 or M.s > 0.5 or M.se > 0.5) {
90
91 if (i < 0 or j + 1 > My - 1 or i + 1 > Mx - 1 or j - 1 < 0) {
92 h_y(i, j, 0) = 0.0;
93 } else {
94 h_y(i, j, 0) = m_h_y_no_model(i, j, 0);
95 }
96
97 }
98
99 // y-component, j-offset
100 if (M.ij > 0.5 or M.n > 0.5) {
101
102 if (j < 0 or j + 1 > My - 1) {
103 h_y(i, j, 1) = 0.0;
104 } else {
105 h_y(i, j, 1) = m_h_y_no_model(i, j, 1);
106 }
107
108 }
109 } // end of the loop over grid points
110 }
111
112 } // end of namespace stressbalance
113
114 } // end of namespace pism