tPico.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
---
tPico.cc (31360B)
---
1 // Copyright (C) 2012-2019 Constantine Khrulev, Ricarda Winkelmann, Ronja Reese, Torsten
2 // Albrecht, and Matthias Mengel
3 //
4 // This file is part of PISM.
5 //
6 // PISM is free software; you can redistribute it and/or modify it under the
7 // terms of the GNU General Public License as published by the Free Software
8 // Foundation; either version 2 of the License, or (at your option) any later
9 // version.
10 //
11 // PISM is distributed in the hope that it will be useful, but WITHOUT ANY
12 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
13 // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
14 // details.
15 //
16 // You should have received a copy of the GNU General Public License
17 // along with PISM; if not, write to the Free Software
18 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
19 //
20 // Please cite this model as:
21 // 1.
22 // Antarctic sub-shelf melt rates via PICO
23 // R. Reese, T. Albrecht, M. Mengel, X. Asay-Davis and R. Winkelmann
24 // The Cryosphere, 12, 1969-1985, (2018)
25 // DOI: 10.5194/tc-12-1969-2018
26 //
27 // 2.
28 // A box model of circulation and melting in ice shelf caverns
29 // D. Olbers & H. Hellmer
30 // Ocean Dynamics (2010), Volume 60, Issue 1, pp 141–153
31 // DOI: 10.1007/s10236-009-0252-z
32
33 #include <gsl/gsl_math.h> // GSL_NAN
34
35 #include "pism/util/ConfigInterface.hh"
36 #include "pism/util/IceGrid.hh"
37 #include "pism/util/Mask.hh"
38 #include "pism/util/Vars.hh"
39 #include "pism/util/iceModelVec.hh"
40 #include "pism/util/Time.hh"
41 #include "pism/geometry/Geometry.hh"
42
43 #include "pism/coupler/util/options.hh"
44
45 #include "Pico.hh"
46 #include "PicoGeometry.hh"
47 #include "PicoPhysics.hh"
48
49 namespace pism {
50 namespace ocean {
51
52 Pico::Pico(IceGrid::ConstPtr g)
53 : CompleteOceanModel(g, std::shared_ptr<OceanModel>()),
54 m_Soc(m_grid, "pico_salinity", WITHOUT_GHOSTS),
55 m_Soc_box0(m_grid, "pico_salinity_box0", WITHOUT_GHOSTS),
56 m_Toc(m_grid, "pico_temperature", WITHOUT_GHOSTS),
57 m_Toc_box0(m_grid, "pico_temperature_box0", WITHOUT_GHOSTS),
58 m_T_star(m_grid, "pico_T_star", WITHOUT_GHOSTS),
59 m_overturning(m_grid, "pico_overturning", WITHOUT_GHOSTS),
60 m_basal_melt_rate(m_grid, "pico_basal_melt_rate", WITH_GHOSTS),
61 m_basin_mask(m_grid, "basins", WITH_GHOSTS),
62 m_geometry(new PicoGeometry(g)) {
63
64 ForcingOptions opt(*m_grid->ctx(), "ocean.pico");
65
66 {
67 unsigned int buffer_size = m_config->get_number("input.forcing.buffer_size");
68 unsigned int evaluations_per_year = m_config->get_number("input.forcing.evaluations_per_year");
69 bool periodic = opt.period > 0;
70
71 File file(m_grid->com, opt.filename, PISM_NETCDF3, PISM_READONLY);
72
73 m_theta_ocean = IceModelVec2T::ForcingField(m_grid,
74 file,
75 "theta_ocean",
76 "", // no standard name
77 buffer_size,
78 evaluations_per_year,
79 periodic,
80 LINEAR);
81
82 m_salinity_ocean = IceModelVec2T::ForcingField(m_grid,
83 file,
84 "salinity_ocean",
85 "", // no standard name
86 buffer_size,
87 evaluations_per_year,
88 periodic,
89 LINEAR);
90 }
91
92 m_theta_ocean->set_attrs("climate_forcing",
93 "potential temperature of the adjacent ocean",
94 "Kelvin", "Kelvin", "", 0);
95
96 m_salinity_ocean->set_attrs("climate_forcing",
97 "salinity of the adjacent ocean",
98 "g/kg", "g/kg", "", 0);
99
100 m_basin_mask.set_attrs("climate_forcing", "mask determines basins for PICO",
101 "", "", "", 0);
102
103 // computed salinity in ocean boxes
104 m_Soc.set_attrs("model_state", "ocean salinity field",
105 "g/kg", "g/kg", "ocean salinity field", 0);
106 m_Soc.metadata().set_number("_FillValue", 0.0);
107
108 // salinity input for box 1
109 m_Soc_box0.set_attrs("model_state", "ocean base salinity field",
110 "g/kg", "g/kg", "", 0);
111 m_Soc_box0.metadata().set_number("_FillValue", 0.0);
112
113 // computed temperature in ocean boxes
114 m_Toc.set_attrs("model_state", "ocean temperature field",
115 "K", "K", "", 0);
116 m_Toc.metadata().set_number("_FillValue", 0.0);
117
118 // temperature input for box 1
119 m_Toc_box0.set_attrs("model_state", "ocean base temperature",
120 "K", "K", "", 0);
121 m_Toc_box0.metadata().set_number("_FillValue", 0.0);
122
123 m_T_star.set_attrs("model_state", "T_star field",
124 "degree C", "degree C", "", 0);
125 m_T_star.metadata().set_number("_FillValue", 0.0);
126
127 m_overturning.set_attrs("model_state", "cavity overturning",
128 "m^3 s-1", "m^3 s-1", "", 0);
129 m_overturning.metadata().set_number("_FillValue", 0.0);
130
131 m_basal_melt_rate.set_attrs("model_state", "PICO sub-shelf melt rate",
132 "m s-1", "m year-1", "", 0);
133 m_basal_melt_rate.metadata().set_number("_FillValue", 0.0);
134
135 m_shelf_base_temperature->metadata().set_number("_FillValue", 0.0);
136
137 m_n_basins = 0;
138
139 m_n_boxes = m_config->get_number("ocean.pico.number_of_boxes");
140 }
141
142
143 Pico::~Pico() {
144 // empty
145 }
146
147
148 void Pico::init_impl(const Geometry &geometry) {
149 (void) geometry;
150
151 m_log->message(2, "* Initializing the Potsdam Ice-shelf Cavity mOdel for the ocean ...\n");
152
153 ForcingOptions opt(*m_grid->ctx(), "ocean.pico");
154
155 m_theta_ocean->init(opt.filename, opt.period, opt.reference_time);
156 m_salinity_ocean->init(opt.filename, opt.period, opt.reference_time);
157
158 m_basin_mask.regrid(opt.filename, CRITICAL);
159
160 // FIXME: m_n_basins is a misnomer
161 m_n_basins = m_basin_mask.max() + 1;
162
163 m_log->message(4, "PICO basin min=%f,max=%f\n", m_basin_mask.min(), m_basin_mask.max());
164
165 PicoPhysics physics(*m_config);
166
167 m_log->message(2, " -Using %d drainage basins and values: \n"
168 " gamma_T= %.2e, overturning_coeff = %.2e... \n",
169 m_n_basins, physics.gamma_T(), physics.overturning_coeff());
170
171 m_log->message(2, " -Depth of continental shelf for computation of temperature and salinity input\n"
172 " is set for whole domain to continental_shelf_depth=%.0f meter\n",
173 physics.continental_shelf_depth());
174
175 // read time-independent data right away:
176 if (m_theta_ocean->n_records() == 1 and m_salinity_ocean->n_records() == 1) {
177 m_theta_ocean->update(m_grid->ctx()->time()->current(), 0.0);
178 m_salinity_ocean->update(m_grid->ctx()->time()->current(), 0.0);
179 }
180 }
181
182 void Pico::define_model_state_impl(const File &output) const {
183
184 m_basin_mask.define(output);
185 m_Soc_box0.define(output);
186 m_Toc_box0.define(output);
187 m_overturning.define(output);
188
189 OceanModel::define_model_state_impl(output);
190 }
191
192 void Pico::write_model_state_impl(const File &output) const {
193
194 m_basin_mask.write(output);
195 m_Soc_box0.write(output);
196 m_Toc_box0.write(output);
197 m_overturning.write(output);
198
199 OceanModel::define_model_state_impl(output);
200 }
201
202 void Pico::update_impl(const Geometry &geometry, double t, double dt) {
203
204 m_theta_ocean->update(t, dt);
205 m_salinity_ocean->update(t, dt);
206
207 m_theta_ocean->average(t, dt);
208 m_salinity_ocean->average(t, dt);
209
210 // set values that will be used outside of floating ice areas
211 {
212 double T_fill_value = m_config->get_number("constants.fresh_water.melting_point_temperature");
213 double Toc_fill_value = m_Toc.metadata().get_number("_FillValue");
214 double Soc_fill_value = m_Soc.metadata().get_number("_FillValue");
215 double M_fill_value = m_basal_melt_rate.metadata().get_number("_FillValue");
216 double O_fill_value = m_overturning.metadata().get_number("_FillValue");
217
218 m_shelf_base_temperature->set(T_fill_value);
219 m_basal_melt_rate.set(M_fill_value);
220 m_Toc.set(Toc_fill_value);
221 m_Soc.set(Soc_fill_value);
222 m_overturning.set(O_fill_value);
223 m_T_star.set(Toc_fill_value);
224 }
225
226 PicoPhysics physics(*m_config);
227
228 const IceModelVec2S &ice_thickness = geometry.ice_thickness;
229 const IceModelVec2CellType &cell_type = geometry.cell_type;
230 const IceModelVec2S &bed_elevation = geometry.bed_elevation;
231
232 // Geometric part of PICO
233 m_geometry->update(bed_elevation, cell_type);
234
235 // FIXME: m_n_shelves is not really the number of shelves.
236 m_n_shelves = m_geometry->ice_shelf_mask().max() + 1;
237
238 // Physical part of PICO
239 {
240
241 // prepare ocean input temperature and salinity
242 {
243 std::vector<double> basin_temperature(m_n_basins);
244 std::vector<double> basin_salinity(m_n_basins);
245
246 compute_ocean_input_per_basin(physics, m_basin_mask, m_geometry->continental_shelf_mask(), *m_salinity_ocean,
247 *m_theta_ocean, basin_temperature, basin_salinity); // per basin
248
249 set_ocean_input_fields(physics, ice_thickness, cell_type, m_basin_mask, m_geometry->ice_shelf_mask(),
250 basin_temperature, basin_salinity, m_Toc_box0, m_Soc_box0); // per shelf
251 }
252
253 // Use the Beckmann-Goosse parameterization to set reasonable values throughout the
254 // domain.
255 beckmann_goosse(physics,
256 ice_thickness, // input
257 m_geometry->ice_shelf_mask(), // input
258 cell_type, // input
259 m_Toc_box0, // input
260 m_Soc_box0, // input
261 m_basal_melt_rate,
262 *m_shelf_base_temperature,
263 m_Toc,
264 m_Soc);
265
266 // In ice shelves, replace Beckmann-Goosse values using the Olbers and Hellmer model.
267 process_box1(physics,
268 ice_thickness, // input
269 m_geometry->ice_shelf_mask(), // input
270 m_geometry->box_mask(), // input
271 m_Toc_box0, // input
272 m_Soc_box0, // input
273 m_basal_melt_rate,
274 *m_shelf_base_temperature,
275 m_T_star,
276 m_Toc,
277 m_Soc,
278 m_overturning);
279
280 process_other_boxes(physics,
281 ice_thickness, // input
282 m_geometry->ice_shelf_mask(), // input
283 m_geometry->box_mask(), // input
284 m_basal_melt_rate,
285 *m_shelf_base_temperature,
286 m_T_star,
287 m_Toc,
288 m_Soc);
289 }
290
291 extend_basal_melt_rates(cell_type,m_basal_melt_rate);
292
293 m_shelf_base_mass_flux->copy_from(m_basal_melt_rate);
294 m_shelf_base_mass_flux->scale(physics.ice_density());
295
296 m_melange_back_pressure_fraction->set(m_config->get_number("ocean.melange_back_pressure_fraction"));
297 }
298
299
300 MaxTimestep Pico::max_timestep_impl(double t) const {
301 (void) t;
302
303 return MaxTimestep("ocean pico");
304 }
305
306
307 //! Compute temperature and salinity input from ocean data by averaging.
308
309 //! We average the ocean data over the continental shelf reagion for each basin.
310 //! We use dummy ocean data if no such average can be calculated.
311
312
313 void Pico::compute_ocean_input_per_basin(const PicoPhysics &physics, const IceModelVec2Int &basin_mask,
314 const IceModelVec2Int &continental_shelf_mask,
315 const IceModelVec2S &salinity_ocean, const IceModelVec2S &theta_ocean,
316 std::vector<double> &temperature, std::vector<double> &salinity) {
317
318 std::vector<int> count(m_n_basins, 0);
319
320 temperature.resize(m_n_basins);
321 salinity.resize(m_n_basins);
322 for (int basin_id = 0; basin_id < m_n_basins; basin_id++) {
323 temperature[basin_id] = 0.0;
324 salinity[basin_id] = 0.0;
325 }
326
327 IceModelVec::AccessList list{ &theta_ocean, &salinity_ocean, &basin_mask, &continental_shelf_mask };
328
329 // compute the sum for each basin for region that intersects with the continental shelf
330 // area and is not covered by an ice shelf. (continental shelf mask excludes ice shelf
331 // areas)
332 for (Points p(*m_grid); p; p.next()) {
333 const int i = p.i(), j = p.j();
334
335 if (continental_shelf_mask.as_int(i, j) == 2) {
336 int basin_id = basin_mask.as_int(i, j);
337
338 count[basin_id] += 1;
339 salinity[basin_id] += salinity_ocean(i, j);
340 temperature[basin_id] += theta_ocean(i, j);
341 }
342 }
343
344 // Divide by number of grid cells if more than zero cells belong to the basin. if no
345 // ocean_contshelf_mask values intersect with the basin, count is zero. In such case,
346 // use dummy temperature and salinity. This could happen, for example, if the ice shelf
347 // front advances beyond the continental shelf break.
348 for (int basin_id = 0; basin_id < m_n_basins; basin_id++) {
349
350 count[basin_id] = GlobalSum(m_grid->com, count[basin_id]);
351 salinity[basin_id] = GlobalSum(m_grid->com, salinity[basin_id]);
352 temperature[basin_id] = GlobalSum(m_grid->com, temperature[basin_id]);
353
354 // if basin is not dummy basin 0 or there are no ocean cells in this basin to take the mean over.
355 if (basin_id > 0 && count[basin_id] == 0) {
356 m_log->message(2, "PICO ocean WARNING: basin %d contains no cells with ocean data on continental shelf\n"
357 "(no values with ocean_contshelf_mask=2).\n"
358 "No mean salinity or temperature values are computed, instead using\n"
359 "the standard values T_dummy =%.3f, S_dummy=%.3f.\n"
360 "This might bias your basal melt rates, check your input data carefully.\n",
361 basin_id, physics.T_dummy(), physics.S_dummy());
362
363 temperature[basin_id] = physics.T_dummy();
364 salinity[basin_id] = physics.S_dummy();
365
366 } else {
367
368 salinity[basin_id] /= count[basin_id];
369 temperature[basin_id] /= count[basin_id];
370
371 m_log->message(5, " %d: temp =%.3f, salinity=%.3f\n", basin_id, temperature[basin_id], salinity[basin_id]);
372 }
373 }
374 }
375
376 //! Set ocean ocean input from box 0 as boundary condition for box 1.
377
378 //! Set ocean temperature and salinity (Toc_box0, Soc_box0)
379 //! from box 0 (in front of the ice shelf) as inputs for
380 //! box 1, which is the ocean box adjacent to the grounding line.
381 //!
382 //! We enforce that Toc_box0 is always at least the local pressure melting point.
383 void Pico::set_ocean_input_fields(const PicoPhysics &physics, const IceModelVec2S &ice_thickness,
384 const IceModelVec2CellType &mask, const IceModelVec2Int &basin_mask,
385 const IceModelVec2Int &shelf_mask, const std::vector<double> basin_temperature,
386 const std::vector<double> basin_salinity, IceModelVec2S &Toc_box0,
387 IceModelVec2S &Soc_box0) {
388
389 IceModelVec::AccessList list{ &ice_thickness, &basin_mask, &Soc_box0, &Toc_box0, &mask, &shelf_mask };
390
391 std::vector<std::vector<int> > n_shelf_cells_per_basin(m_n_shelves, std::vector<int>(m_n_basins, 0));
392 std::vector<int> n_shelf_cells(m_n_shelves, 0);
393
394 // 1) count the number of cells in each shelf
395 // 2) count the number of cells in the intersection of each shelf with all the basins
396 {
397 for (Points p(*m_grid); p; p.next()) {
398 const int i = p.i(), j = p.j();
399 int s = shelf_mask.as_int(i, j);
400 int b = basin_mask.as_int(i, j);
401 n_shelf_cells_per_basin[s][b]++;
402 n_shelf_cells[s]++;
403 }
404
405 for (int s = 0; s < m_n_shelves; s++) {
406 n_shelf_cells[s] = GlobalSum(m_grid->com, n_shelf_cells[s]);
407 for (int b = 0; b < m_n_basins; b++) {
408 n_shelf_cells_per_basin[s][b] = GlobalSum(m_grid->com, n_shelf_cells_per_basin[s][b]);
409 }
410 }
411 }
412
413 // now set potential temperature and salinity box 0:
414
415 int low_temperature_counter = 0;
416 for (Points p(*m_grid); p; p.next()) {
417 const int i = p.i(), j = p.j();
418
419 // make sure all temperatures are zero at the beginning of each time step
420 Toc_box0(i, j) = 0.0; // in K
421 Soc_box0(i, j) = 0.0; // in psu
422
423 int s = shelf_mask.as_int(i, j);
424
425 if (mask.as_int(i, j) == MASK_FLOATING and s > 0) {
426 // note: shelf_mask = 0 in lakes
427
428 // weighted input depending on the number of shelf cells in each basin
429 for (int b = 1; b < m_n_basins; b++) { //Note: b=0 yields nan
430 Toc_box0(i, j) += basin_temperature[b] * n_shelf_cells_per_basin[s][b] / (double)n_shelf_cells[s];
431 Soc_box0(i, j) += basin_salinity[b] * n_shelf_cells_per_basin[s][b] / (double)n_shelf_cells[s];
432 }
433
434 double theta_pm = physics.theta_pm(Soc_box0(i, j), physics.pressure(ice_thickness(i, j)));
435
436 // temperature input for grounding line box should not be below pressure melting point
437 if (Toc_box0(i, j) < theta_pm) {
438 // Setting Toc_box0 a little higher than theta_pm ensures that later equations are well solvable.
439 Toc_box0(i, j) = theta_pm + 0.001;
440 low_temperature_counter += 1;
441 }
442 }
443 }
444
445 low_temperature_counter = GlobalSum(m_grid->com, low_temperature_counter);
446 if (low_temperature_counter > 0) {
447 m_log->message(2, "PICO ocean warning: temperature has been below pressure melting temperature in %d cases,\n"
448 "setting it to pressure melting temperature\n",
449 low_temperature_counter);
450 }
451 }
452
453 /*!
454 * Use the simpler parameterization due to [@ref BeckmannGoosse2003] to set default
455 * sub-shelf temperature and melt rate values.
456 *
457 * At grid points containing floating ice not connected to the ocean, set the basal melt
458 * rate to zero and set basal temperature to the pressure melting point.
459 */
460 void Pico::beckmann_goosse(const PicoPhysics &physics,
461 const IceModelVec2S &ice_thickness,
462 const IceModelVec2Int &shelf_mask,
463 const IceModelVec2CellType &cell_type,
464 const IceModelVec2S &Toc_box0,
465 const IceModelVec2S &Soc_box0,
466 IceModelVec2S &basal_melt_rate,
467 IceModelVec2S &basal_temperature,
468 IceModelVec2S &Toc,
469 IceModelVec2S &Soc) {
470
471 const double T0 = m_config->get_number("constants.fresh_water.melting_point_temperature"),
472 beta_CC = m_config->get_number("constants.ice.beta_Clausius_Clapeyron"),
473 g = m_config->get_number("constants.standard_gravity"),
474 ice_density = m_config->get_number("constants.ice.density");
475
476 IceModelVec::AccessList list{ &ice_thickness, &cell_type, &shelf_mask, &Toc_box0, &Soc_box0,
477 &Toc, &Soc, &basal_melt_rate, &basal_temperature };
478
479 for (Points p(*m_grid); p; p.next()) {
480 const int i = p.i(), j = p.j();
481
482 if (cell_type.floating_ice(i, j)) {
483 if (shelf_mask.as_int(i, j) > 0) {
484 double pressure = physics.pressure(ice_thickness(i, j));
485
486 basal_melt_rate(i, j) =
487 physics.melt_rate_beckmann_goosse(physics.theta_pm(Soc_box0(i, j), pressure), Toc_box0(i, j));
488 basal_temperature(i, j) = physics.T_pm(Soc_box0(i, j), pressure);
489
490 // diagnostic outputs
491 Toc(i, j) = Toc_box0(i, j); // in Kelvin
492 Soc(i, j) = Soc_box0(i, j); // in psu
493 } else {
494 // Floating ice cells not connected to the ocean.
495 const double pressure = ice_density * g * ice_thickness(i, j); // FIXME issue #15
496
497 basal_temperature(i, j) = T0 - beta_CC * pressure;
498 basal_melt_rate(i, j) = 0.0;
499 }
500 }
501 }
502 }
503
504
505 void Pico::process_box1(const PicoPhysics &physics,
506 const IceModelVec2S &ice_thickness,
507 const IceModelVec2Int &shelf_mask,
508 const IceModelVec2Int &box_mask,
509 const IceModelVec2S &Toc_box0,
510 const IceModelVec2S &Soc_box0,
511 IceModelVec2S &basal_melt_rate,
512 IceModelVec2S &basal_temperature,
513 IceModelVec2S &T_star,
514 IceModelVec2S &Toc,
515 IceModelVec2S &Soc,
516 IceModelVec2S &overturning) {
517
518 std::vector<double> box1_area(m_n_shelves);
519
520 compute_box_area(1, shelf_mask, box_mask, box1_area);
521
522 IceModelVec::AccessList list{ &ice_thickness, &shelf_mask, &box_mask, &T_star, &Toc_box0, &Toc,
523 &Soc_box0, &Soc, &overturning, &basal_melt_rate, &basal_temperature };
524
525 int n_Toc_failures = 0;
526
527 // basal melt rate, ambient temperature and salinity and overturning calculation
528 // for each box1 grid cell.
529 for (Points p(*m_grid); p; p.next()) {
530 const int i = p.i(), j = p.j();
531
532 int shelf_id = shelf_mask.as_int(i, j);
533
534 if (box_mask.as_int(i, j) == 1 and shelf_id > 0) {
535
536 const double pressure = physics.pressure(ice_thickness(i, j));
537
538 T_star(i, j) = physics.T_star(Soc_box0(i, j), Toc_box0(i, j), pressure);
539
540 auto Toc_box1 = physics.Toc_box1(box1_area[shelf_id], T_star(i, j), Soc_box0(i, j), Toc_box0(i, j));
541
542 // This can only happen if T_star > 0.25*p_coeff, in particular T_star > 0
543 // which can only happen for values of Toc_box0 close to the local pressure melting point
544 if (Toc_box1.failed) {
545 m_log->message(5, "PICO ocean WARNING: negative square root argument at %d, %d\n"
546 "probably because of positive T_star=%f \n"
547 "Not aborting, but setting square root to 0... \n",
548 i, j, T_star(i, j));
549
550 n_Toc_failures += 1;
551 }
552
553 Toc(i, j) = Toc_box1.value;
554 Soc(i, j) = physics.Soc_box1(Toc_box0(i, j), Soc_box0(i, j), Toc(i, j)); // in psu
555
556 overturning(i, j) = physics.overturning(Soc_box0(i, j), Soc(i, j), Toc_box0(i, j), Toc(i, j));
557
558 // main outputs
559 basal_melt_rate(i, j) = physics.melt_rate(physics.theta_pm(Soc(i, j), pressure), Toc(i, j));
560 basal_temperature(i, j) = physics.T_pm(Soc(i, j), pressure);
561 }
562 }
563
564 n_Toc_failures = GlobalSum(m_grid->com, n_Toc_failures);
565 if (n_Toc_failures > 0) {
566 m_log->message(2, "PICO ocean warning: square-root argument for temperature calculation "
567 "has been negative in %d cases!\n",
568 n_Toc_failures);
569 }
570 }
571
572 void Pico::process_other_boxes(const PicoPhysics &physics,
573 const IceModelVec2S &ice_thickness,
574 const IceModelVec2Int &shelf_mask,
575 const IceModelVec2Int &box_mask,
576 IceModelVec2S &basal_melt_rate,
577 IceModelVec2S &basal_temperature,
578 IceModelVec2S &T_star,
579 IceModelVec2S &Toc,
580 IceModelVec2S &Soc) {
581
582 std::vector<double> overturning(m_n_shelves, 0.0);
583 std::vector<double> salinity(m_n_shelves, 0.0);
584 std::vector<double> temperature(m_n_shelves, 0.0);
585
586 // get average overturning from box 1 that is used as input later
587 compute_box_average(1, m_overturning, shelf_mask, box_mask, overturning);
588
589 std::vector<bool> use_beckmann_goosse(m_n_shelves);
590
591 IceModelVec::AccessList list{ &ice_thickness, &shelf_mask, &box_mask, &T_star, &Toc,
592 &Soc, &basal_melt_rate, &basal_temperature };
593
594 // Iterate over all boxes i for i > 1
595 for (int box = 2; box <= m_n_boxes; ++box) {
596
597 compute_box_average(box - 1, Toc, shelf_mask, box_mask, temperature);
598 compute_box_average(box - 1, Soc, shelf_mask, box_mask, salinity);
599
600 // find all the shelves where we should fall back to the Beckmann-Goosse
601 // parameterization
602 for (int s = 1; s < m_n_shelves; ++s) {
603 if (salinity[s] == 0.0 or temperature[s] == 0.0 or overturning[s] == 0.0) {
604 use_beckmann_goosse[s] = true;
605 } else {
606 use_beckmann_goosse[s] = false;
607 }
608 }
609
610 std::vector<double> box_area;
611 compute_box_area(box, shelf_mask, box_mask, box_area);
612
613 int n_beckmann_goosse_cells = 0;
614
615 for (Points p(*m_grid); p; p.next()) {
616 const int i = p.i(), j = p.j();
617
618 int shelf_id = shelf_mask.as_int(i, j);
619
620 if (box_mask.as_int(i, j) == box and shelf_id > 0) {
621
622 if (use_beckmann_goosse[shelf_id]) {
623 n_beckmann_goosse_cells += 1;
624 continue;
625 }
626
627 // Get the input from previous box
628 const double
629 S_previous = salinity[shelf_id],
630 T_previous = temperature[shelf_id],
631 overturning_box1 = overturning[shelf_id];
632
633 {
634 double pressure = physics.pressure(ice_thickness(i, j));
635
636 // diagnostic outputs
637 T_star(i, j) = physics.T_star(S_previous, T_previous, pressure);
638 Toc(i, j) = physics.Toc(box_area[shelf_id], T_previous, T_star(i, j), overturning_box1, S_previous);
639 Soc(i, j) = physics.Soc(S_previous, T_previous, Toc(i, j));
640
641 // main outputs: basal melt rate and temperature
642 basal_melt_rate(i, j) = physics.melt_rate(physics.theta_pm(Soc(i, j), pressure), Toc(i, j));
643 basal_temperature(i, j) = physics.T_pm(Soc(i, j), pressure);
644 }
645 }
646 } // loop over grid points
647
648 n_beckmann_goosse_cells = GlobalSum(m_grid->com, n_beckmann_goosse_cells);
649 if (n_beckmann_goosse_cells > 0) {
650 m_log->message(2, "PICO ocean WARNING: box %d, no boundary data from previous box in %d case(s)!\n"
651 "switching to Beckmann Goosse (2003) meltrate calculation\n",
652 box, n_beckmann_goosse_cells);
653 }
654
655 } // loop over boxes
656 }
657
658 /*!
659 * Extend basal melt rates to grounded and ocean neighbors for consitency with subgl_melt.
660 * Note that melt rates are then simply interpolated into partially floating cells, they
661 * are not included in the calculations of PICO.
662 */
663 void Pico::extend_basal_melt_rates(const IceModelVec2CellType &cell_type, IceModelVec2S &basal_melt_rate) {
664
665 IceModelVec::AccessList list{&cell_type, &basal_melt_rate};
666
667 for (Points p(*m_grid); p; p.next()) {
668
669 const int i = p.i(), j = p.j();
670
671 auto M = cell_type.int_box(i, j);
672
673 bool potential_partially_filled_cell =
674 ((M.ij == MASK_GROUNDED or M.ij == MASK_ICE_FREE_OCEAN) and
675 (M.w == MASK_FLOATING or M.e == MASK_FLOATING or M.s == MASK_FLOATING or M.n == MASK_FLOATING or
676 M.sw == MASK_FLOATING or M.nw == MASK_FLOATING or M.se == MASK_FLOATING or M.ne == MASK_FLOATING) );
677
678 if (potential_partially_filled_cell) {
679 auto BMR = basal_melt_rate.box(i, j);
680
681 int N = 0;
682 double melt_sum = 0.0;
683
684 melt_sum += M.nw == MASK_FLOATING ? (++N, BMR.nw) : 0.0;
685 melt_sum += M.n == MASK_FLOATING ? (++N, BMR.n) : 0.0;
686 melt_sum += M.ne == MASK_FLOATING ? (++N, BMR.ne) : 0.0;
687 melt_sum += M.e == MASK_FLOATING ? (++N, BMR.e) : 0.0;
688 melt_sum += M.se == MASK_FLOATING ? (++N, BMR.se) : 0.0;
689 melt_sum += M.s == MASK_FLOATING ? (++N, BMR.s) : 0.0;
690 melt_sum += M.sw == MASK_FLOATING ? (++N, BMR.sw) : 0.0;
691 melt_sum += M.w == MASK_FLOATING ? (++N, BMR.w) : 0.0;
692
693 if (N != 0) { // If there are floating neigbors, return average melt rates
694 basal_melt_rate(i, j) = melt_sum / N;
695 }
696 }
697 } // end of the loop over grid points
698 }
699
700
701
702 // Write diagnostic variables to extra files if requested
703 DiagnosticList Pico::diagnostics_impl() const {
704
705 DiagnosticList result = {
706 { "basins", Diagnostic::wrap(m_basin_mask) },
707 { "pico_overturning", Diagnostic::wrap(m_overturning) },
708 { "pico_salinity_box0", Diagnostic::wrap(m_Soc_box0) },
709 { "pico_temperature_box0", Diagnostic::wrap(m_Toc_box0) },
710 { "pico_box_mask", Diagnostic::wrap(m_geometry->box_mask()) },
711 { "pico_shelf_mask", Diagnostic::wrap(m_geometry->ice_shelf_mask()) },
712 { "pico_ice_rise_mask", Diagnostic::wrap(m_geometry->ice_rise_mask()) },
713 { "pico_basal_melt_rate", Diagnostic::wrap(m_basal_melt_rate) },
714 { "pico_contshelf_mask", Diagnostic::wrap(m_geometry->continental_shelf_mask()) },
715 { "pico_salinity", Diagnostic::wrap(m_Soc) },
716 { "pico_temperature", Diagnostic::wrap(m_Toc) },
717 { "pico_T_star", Diagnostic::wrap(m_T_star) },
718 { "pico_basal_temperature", Diagnostic::wrap(*m_shelf_base_temperature) },
719 };
720
721 return combine(result, OceanModel::diagnostics_impl());
722 }
723
724 /*!
725 * For each shelf, compute average of a given field over the box with id `box_id`.
726 *
727 * This method is used to get inputs from a previous box for the next one.
728 */
729 void Pico::compute_box_average(int box_id,
730 const IceModelVec2S &field,
731 const IceModelVec2Int &shelf_mask,
732 const IceModelVec2Int &box_mask,
733 std::vector<double> &result) {
734
735 IceModelVec::AccessList list{ &field, &shelf_mask, &box_mask };
736
737 std::vector<int> n_cells_per_box(m_n_shelves, 0);
738
739 // fill results with zeros
740 result.resize(m_n_shelves);
741 for (int s = 0; s < m_n_shelves; ++s) {
742 result[s] = 0.0;
743 }
744
745 // compute the sum of field in each shelf's box box_id
746 for (Points p(*m_grid); p; p.next()) {
747 const int i = p.i(), j = p.j();
748
749 int shelf_id = shelf_mask.as_int(i, j);
750
751 if (box_mask.as_int(i, j) == box_id) {
752 n_cells_per_box[shelf_id] += 1;
753 result[shelf_id] += field(i, j);
754 }
755 }
756
757 // compute the global sum and average
758 for (int s = 0; s < m_n_shelves; ++s) {
759 auto n_cells = GlobalSum(m_grid->com, n_cells_per_box[s]);
760
761 result[s] = GlobalSum(m_grid->com, result[s]);
762
763 if (n_cells > 0) {
764 result[s] /= (double)n_cells;
765 }
766 }
767 }
768
769 /*!
770 * For all shelves compute areas of boxes with id `box_id`.
771 *
772 * @param[in] box_mask box index mask
773 * @param[out] result resulting box areas.
774 *
775 * Note: shelf and box indexes start from 1.
776 */
777 void Pico::compute_box_area(int box_id,
778 const IceModelVec2Int &shelf_mask,
779 const IceModelVec2Int &box_mask,
780 std::vector<double> &result) {
781 result.resize(m_n_shelves);
782
783 IceModelVec::AccessList list{ &shelf_mask, &box_mask };
784
785 auto cell_area = m_grid->cell_area();
786
787 for (Points p(*m_grid); p; p.next()) {
788 const int i = p.i(), j = p.j();
789
790 int shelf_id = shelf_mask.as_int(i, j);
791
792 if (shelf_id > 0 and box_mask.as_int(i, j) == box_id) {
793 result[shelf_id] += cell_area;
794 }
795 }
796
797 // compute global sums
798 for (int s = 1; s < m_n_shelves; ++s) {
799 result[s] = GlobalSum(m_grid->com, result[s]);
800 }
801 }
802
803 } // end of namespace ocean
804 } // end of namespace pism