PISM, A Parallel Ice Sheet Model  stable v2.0.6 committed by Constantine Khrulev on 2023-01-23 15:14:38 -0900
Pico.cc
Go to the documentation of this file.
1 // Copyright (C) 2012-2019, 2021 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/coupler/util/options.hh"
36 #include "pism/geometry/Geometry.hh"
37 #include "pism/util/ConfigInterface.hh"
38 #include "pism/util/IceGrid.hh"
39 #include "pism/util/Mask.hh"
40 #include "pism/util/Time.hh"
41 #include "pism/util/Vars.hh"
42 #include "pism/util/iceModelVec.hh"
43 
44 #include "Pico.hh"
45 #include "PicoGeometry.hh"
46 #include "PicoPhysics.hh"
47 
48 namespace pism {
49 namespace ocean {
50 
52  : CompleteOceanModel(grid, std::shared_ptr<OceanModel>()),
53  m_Soc(m_grid, "pico_salinity", WITHOUT_GHOSTS),
54  m_Soc_box0(m_grid, "pico_salinity_box0", WITHOUT_GHOSTS),
55  m_Toc(m_grid, "pico_temperature", WITHOUT_GHOSTS),
56  m_Toc_box0(m_grid, "pico_temperature_box0", WITHOUT_GHOSTS),
57  m_T_star(m_grid, "pico_T_star", WITHOUT_GHOSTS),
58  m_overturning(m_grid, "pico_overturning", WITHOUT_GHOSTS),
59  m_basal_melt_rate(m_grid, "pico_basal_melt_rate", WITH_GHOSTS),
60  m_geometry(grid),
61  m_n_basins(0),
62  m_n_boxes(0),
63  m_n_shelves(0) {
64 
65  ForcingOptions opt(*m_grid->ctx(), "ocean.pico");
66 
67  {
68  auto buffer_size = static_cast<int>(m_config->get_number("input.forcing.buffer_size"));
69 
70  File file(m_grid->com, opt.filename, PISM_NETCDF3, PISM_READONLY);
71 
73  file,
74  "theta_ocean",
75  "", // no standard name
76  buffer_size,
77  opt.periodic,
78  LINEAR);
79 
81  file,
82  "salinity_ocean",
83  "", // no standard name
84  buffer_size,
85  opt.periodic,
86  LINEAR);
87  }
88 
89  m_theta_ocean->set_attrs("climate_forcing",
90  "potential temperature of the adjacent ocean",
91  "Kelvin", "Kelvin", "", 0);
92 
93  m_salinity_ocean->set_attrs("climate_forcing",
94  "salinity of the adjacent ocean",
95  "g/kg", "g/kg", "", 0);
96 
97  // computed salinity in ocean boxes
98  m_Soc.set_attrs("model_state", "ocean salinity field",
99  "g/kg", "g/kg", "ocean salinity field", 0);
100  m_Soc.metadata()["_FillValue"] = {0.0};
101 
102  // salinity input for box 1
103  m_Soc_box0.set_attrs("model_state", "ocean base salinity field",
104  "g/kg", "g/kg", "", 0);
105  m_Soc_box0.metadata()["_FillValue"] = {0.0};
106 
107  // computed temperature in ocean boxes
108  m_Toc.set_attrs("model_state", "ocean temperature field",
109  "K", "K", "", 0);
110  m_Toc.metadata()["_FillValue"] = {0.0};
111 
112  // temperature input for box 1
113  m_Toc_box0.set_attrs("model_state", "ocean base temperature",
114  "K", "K", "", 0);
115  m_Toc_box0.metadata()["_FillValue"] = {0.0};
116 
117  m_T_star.set_attrs("model_state", "T_star field",
118  "degree C", "degree C", "", 0);
119  m_T_star.metadata()["_FillValue"] = {0.0};
120 
121  m_overturning.set_attrs("model_state", "cavity overturning",
122  "m^3 s-1", "m^3 s-1", "", 0);
123  m_overturning.metadata()["_FillValue"] = {0.0};
124 
125  m_basal_melt_rate.set_attrs("model_state", "PICO sub-shelf melt rate",
126  "m s-1", "m year-1", "", 0);
127  m_basal_melt_rate.metadata()["_FillValue"] = {0.0};
128 
129  m_shelf_base_temperature->metadata()["_FillValue"] = {0.0};
130 
131  m_n_boxes = static_cast<int>(m_config->get_number("ocean.pico.number_of_boxes"));
132 }
133 
134 void Pico::init_impl(const Geometry &geometry) {
135  (void) geometry;
136 
137  m_log->message(2, "* Initializing the Potsdam Ice-shelf Cavity mOdel for the ocean ...\n");
138 
139  ForcingOptions opt(*m_grid->ctx(), "ocean.pico");
140 
141  m_theta_ocean->init(opt.filename, opt.periodic);
142  m_salinity_ocean->init(opt.filename, opt.periodic);
143 
144  // This initializes the basin_mask
145  m_geometry.init();
146 
147  // FIXME: m_n_basins is a misnomer
148  m_n_basins = static_cast<int>(max(m_geometry.basin_mask())) + 1;
149 
150  m_log->message(4, "PICO basin min=%f, max=%f\n",
153 
154  PicoPhysics physics(*m_config);
155 
156  m_log->message(2, " -Using %d drainage basins and values: \n"
157  " gamma_T= %.2e, overturning_coeff = %.2e... \n",
158  m_n_basins, physics.gamma_T(), physics.overturning_coeff());
159 
160  m_log->message(2, " -Depth of continental shelf for computation of temperature and salinity input\n"
161  " is set for whole domain to continental_shelf_depth=%.0f meter\n",
162  physics.continental_shelf_depth());
163 
164  // read time-independent data right away:
165  if (m_theta_ocean->buffer_size() == 1 and m_salinity_ocean->buffer_size() == 1) {
166  m_theta_ocean->update(m_grid->ctx()->time()->current(), 0.0);
167  m_salinity_ocean->update(m_grid->ctx()->time()->current(), 0.0);
168  }
169 
170  double
171  ice_density = m_config->get_number("constants.ice.density"),
172  water_density = m_config->get_number("constants.sea_water.density"),
173  g = m_config->get_number("constants.standard_gravity");
174 
175  compute_average_water_column_pressure(geometry, ice_density, water_density, g,
177 }
178 
179 void Pico::define_model_state_impl(const File &output) const {
180 
181  m_geometry.basin_mask().define(output);
182  m_Soc_box0.define(output);
183  m_Toc_box0.define(output);
184  m_overturning.define(output);
185 
187 }
188 
189 void Pico::write_model_state_impl(const File &output) const {
190 
191  m_geometry.basin_mask().write(output);
192  m_Soc_box0.write(output);
193  m_Toc_box0.write(output);
194  m_overturning.write(output);
195 
197 }
198 
199 /*!
200 * Extend basal melt rates to grounded and ocean neighbors for consitency with subgl_melt.
201 * Note that melt rates are then simply interpolated into partially floating cells, they
202 * are not included in the calculations of PICO.
203 */
204 static void extend_basal_melt_rates(const IceModelVec2CellType &cell_type,
205  IceModelVec2S &basal_melt_rate) {
206 
207  auto grid = basal_melt_rate.grid();
208 
209  // update ghosts of the basal melt rate so that we can use basal_melt_rate.box(i,j)
210  // below
211  basal_melt_rate.update_ghosts();
212 
213  IceModelVec::AccessList list{&cell_type, &basal_melt_rate};
214 
215  for (Points p(*grid); p; p.next()) {
216 
217  const int i = p.i(), j = p.j();
218 
219  auto M = cell_type.box(i, j);
220 
221  bool potential_partially_filled_cell =
222  ((M.ij == MASK_GROUNDED or M.ij == MASK_ICE_FREE_OCEAN) and
223  (M.w == MASK_FLOATING or M.e == MASK_FLOATING or
224  M.s == MASK_FLOATING or M.n == MASK_FLOATING or
225  M.sw == MASK_FLOATING or M.nw == MASK_FLOATING or
226  M.se == MASK_FLOATING or M.ne == MASK_FLOATING));
227 
228  if (potential_partially_filled_cell) {
229  auto BMR = basal_melt_rate.box(i, j);
230 
231  int N = 0;
232  double melt_sum = 0.0;
233 
234  melt_sum += M.nw == MASK_FLOATING ? (++N, BMR.nw) : 0.0;
235  melt_sum += M.n == MASK_FLOATING ? (++N, BMR.n) : 0.0;
236  melt_sum += M.ne == MASK_FLOATING ? (++N, BMR.ne) : 0.0;
237  melt_sum += M.e == MASK_FLOATING ? (++N, BMR.e) : 0.0;
238  melt_sum += M.se == MASK_FLOATING ? (++N, BMR.se) : 0.0;
239  melt_sum += M.s == MASK_FLOATING ? (++N, BMR.s) : 0.0;
240  melt_sum += M.sw == MASK_FLOATING ? (++N, BMR.sw) : 0.0;
241  melt_sum += M.w == MASK_FLOATING ? (++N, BMR.w) : 0.0;
242 
243  if (N != 0) { // If there are floating neigbors, return average melt rates
244  basal_melt_rate(i, j) = melt_sum / N;
245  }
246  }
247  } // end of the loop over grid points
248 }
249 
250 void Pico::update_impl(const Geometry &geometry, double t, double dt) {
251 
252  m_theta_ocean->update(t, dt);
253  m_salinity_ocean->update(t, dt);
254 
255  m_theta_ocean->average(t, dt);
256  m_salinity_ocean->average(t, dt);
257 
258  // set values that will be used outside of floating ice areas
259  {
260  double T_fill_value = m_config->get_number("constants.fresh_water.melting_point_temperature");
261  double Toc_fill_value = m_Toc.metadata().get_number("_FillValue");
262  double Soc_fill_value = m_Soc.metadata().get_number("_FillValue");
263  double M_fill_value = m_basal_melt_rate.metadata().get_number("_FillValue");
264  double O_fill_value = m_overturning.metadata().get_number("_FillValue");
265 
266  m_shelf_base_temperature->set(T_fill_value);
267  m_basal_melt_rate.set(M_fill_value);
268  m_Toc.set(Toc_fill_value);
269  m_Soc.set(Soc_fill_value);
270  m_overturning.set(O_fill_value);
271  m_T_star.set(Toc_fill_value);
272  }
273 
274  PicoPhysics physics(*m_config);
275 
276  const IceModelVec2S &ice_thickness = geometry.ice_thickness;
277  const IceModelVec2CellType &cell_type = geometry.cell_type;
278  const IceModelVec2S &bed_elevation = geometry.bed_elevation;
279 
280  // Geometric part of PICO
281  m_geometry.update(bed_elevation, cell_type);
282 
283  // FIXME: m_n_shelves is not really the number of shelves.
284  m_n_shelves = static_cast<int>(max(m_geometry.ice_shelf_mask())) + 1;
285 
286  // Physical part of PICO
287  {
288 
289  // prepare ocean input temperature and salinity
290  {
291  std::vector<double> basin_temperature(m_n_basins);
292  std::vector<double> basin_salinity(m_n_basins);
293 
298  *m_theta_ocean,
299  basin_temperature,
300  basin_salinity); // per basin
301 
302  set_ocean_input_fields(physics,
303  ice_thickness,
304  cell_type,
307  basin_temperature,
308  basin_salinity,
309  m_Toc_box0,
310  m_Soc_box0); // per shelf
311  }
312 
313  // Use the Beckmann-Goosse parameterization to set reasonable values throughout the
314  // domain.
315  beckmann_goosse(physics,
316  ice_thickness, // input
317  m_geometry.ice_shelf_mask(), // input
318  cell_type, // input
319  m_Toc_box0, // input
320  m_Soc_box0, // input
323  m_Toc,
324  m_Soc);
325 
326  // In ice shelves, replace Beckmann-Goosse values using the Olbers and Hellmer model.
327  process_box1(physics,
328  ice_thickness, // input
329  m_geometry.ice_shelf_mask(), // input
330  m_geometry.box_mask(), // input
331  m_Toc_box0, // input
332  m_Soc_box0, // input
335  m_T_star,
336  m_Toc,
337  m_Soc,
338  m_overturning);
339 
340  process_other_boxes(physics,
341  ice_thickness, // input
342  m_geometry.ice_shelf_mask(), // input
343  m_geometry.box_mask(), // input
346  m_T_star,
347  m_Toc,
348  m_Soc);
349  }
350 
352 
354  m_shelf_base_mass_flux->scale(physics.ice_density());
355 
356  double
357  ice_density = m_config->get_number("constants.ice.density"),
358  water_density = m_config->get_number("constants.sea_water.density"),
359  g = m_config->get_number("constants.standard_gravity");
360 
361  compute_average_water_column_pressure(geometry, ice_density, water_density, g,
363 }
364 
365 
367  (void) t;
368 
369  return MaxTimestep("ocean pico");
370 }
371 
372 
373 //! Compute temperature and salinity input from ocean data by averaging.
374 
375 //! We average the ocean data over the continental shelf reagion for each basin.
376 //! We use dummy ocean data if no such average can be calculated.
377 
378 
380  const IceModelVec2Int &basin_mask,
381  const IceModelVec2Int &continental_shelf_mask,
382  const IceModelVec2S &salinity_ocean,
383  const IceModelVec2S &theta_ocean,
384  std::vector<double> &temperature,
385  std::vector<double> &salinity) const {
386  std::vector<int> count(m_n_basins, 0);
387  // additional vectors to allreduce efficiently with IntelMPI
388  std::vector<int> countr(m_n_basins, 0);
389  std::vector<double> salinityr(m_n_basins);
390  std::vector<double> temperaturer(m_n_basins);
391 
392  temperature.resize(m_n_basins);
393  salinity.resize(m_n_basins);
394  for (int basin_id = 0; basin_id < m_n_basins; basin_id++) {
395  temperature[basin_id] = 0.0;
396  salinity[basin_id] = 0.0;
397  }
398 
399  IceModelVec::AccessList list{ &theta_ocean, &salinity_ocean, &basin_mask, &continental_shelf_mask };
400 
401  // compute the sum for each basin for region that intersects with the continental shelf
402  // area and is not covered by an ice shelf. (continental shelf mask excludes ice shelf
403  // areas)
404  for (Points p(*m_grid); p; p.next()) {
405  const int i = p.i(), j = p.j();
406 
407  if (continental_shelf_mask.as_int(i, j) == 2) {
408  int basin_id = basin_mask.as_int(i, j);
409 
410  count[basin_id] += 1;
411  salinity[basin_id] += salinity_ocean(i, j);
412  temperature[basin_id] += theta_ocean(i, j);
413  }
414  }
415 
416  // Divide by number of grid cells if more than zero cells belong to the basin. if no
417  // ocean_contshelf_mask values intersect with the basin, count is zero. In such case,
418  // use dummy temperature and salinity. This could happen, for example, if the ice shelf
419  // front advances beyond the continental shelf break.
420  GlobalSum(m_grid->com, count.data(), countr.data(), m_n_basins);
421  GlobalSum(m_grid->com, salinity.data(), salinityr.data(), m_n_basins);
422  GlobalSum(m_grid->com, temperature.data(), temperaturer.data(), m_n_basins);
423 
424  // copy values
425  count = countr;
426  salinity = salinityr;
427  temperature = temperaturer;
428 
429  // "dummy" basin
430  {
431  temperature[0] = physics.T_dummy();
432  salinity[0] = physics.S_dummy();
433  }
434 
435  for (int basin_id = 1; basin_id < m_n_basins; basin_id++) {
436 
437  if (count[basin_id] != 0) {
438  salinity[basin_id] /= count[basin_id];
439  temperature[basin_id] /= count[basin_id];
440 
441  m_log->message(5, " %d: temp =%.3f, salinity=%.3f\n", basin_id, temperature[basin_id], salinity[basin_id]);
442  } else {
443  m_log->message(2, "PICO ocean WARNING: basin %d contains no cells with ocean data on continental shelf\n"
444  "(no values with ocean_contshelf_mask=2).\n"
445  "No mean salinity or temperature values are computed, instead using\n"
446  "the standard values T_dummy =%.3f, S_dummy=%.3f.\n"
447  "This might bias your basal melt rates, check your input data carefully.\n",
448  basin_id, physics.T_dummy(), physics.S_dummy());
449 
450  temperature[basin_id] = physics.T_dummy();
451  salinity[basin_id] = physics.S_dummy();
452  }
453  }
454 }
455 
456 //! Set ocean ocean input from box 0 as boundary condition for box 1.
457 
458 //! Set ocean temperature and salinity (Toc_box0, Soc_box0)
459 //! from box 0 (in front of the ice shelf) as inputs for
460 //! box 1, which is the ocean box adjacent to the grounding line.
461 //!
462 //! We enforce that Toc_box0 is always at least the local pressure melting point.
464  const IceModelVec2S &ice_thickness,
465  const IceModelVec2CellType &mask,
466  const IceModelVec2Int &basin_mask,
467  const IceModelVec2Int &shelf_mask,
468  const std::vector<double> &basin_temperature,
469  const std::vector<double> &basin_salinity,
470  IceModelVec2S &Toc_box0,
471  IceModelVec2S &Soc_box0) const {
472 
473  IceModelVec::AccessList list{ &ice_thickness, &basin_mask, &Soc_box0, &Toc_box0, &mask, &shelf_mask };
474 
475  std::vector<int> n_shelf_cells_per_basin(m_n_shelves * m_n_basins, 0);
476  std::vector<int> n_shelf_cells(m_n_shelves, 0);
477  std::vector<int> cfs_in_basins_per_shelf(m_n_shelves * m_n_basins, 0);
478  // additional vectors to allreduce efficiently with IntelMPI
479  std::vector<int> n_shelf_cells_per_basinr(m_n_shelves * m_n_basins, 0);
480  std::vector<int> n_shelf_cellsr(m_n_shelves, 0);
481  std::vector<int> cfs_in_basins_per_shelfr(m_n_shelves * m_n_basins, 0);
482 
483  // 1) count the number of cells in each shelf
484  // 2) count the number of cells in the intersection of each shelf with all the basins
485  {
486  for (Points p(*m_grid); p; p.next()) {
487  const int i = p.i(), j = p.j();
488  int s = shelf_mask.as_int(i, j);
489  int b = basin_mask.as_int(i, j);
490  n_shelf_cells_per_basin[s*m_n_basins+b]++;
491  n_shelf_cells[s]++;
492 
493  // find all basins b, in which the ice shelf s has a calving front with potential ocean water intrusion
494  if (mask.as_int(i, j) == MASK_FLOATING) {
495  auto M = mask.star(i, j);
496  if (M.n == MASK_ICE_FREE_OCEAN or
497  M.e == MASK_ICE_FREE_OCEAN or
498  M.s == MASK_ICE_FREE_OCEAN or
499  M.w == MASK_ICE_FREE_OCEAN) {
500  if (cfs_in_basins_per_shelf[s * m_n_basins + b] != b) {
501  cfs_in_basins_per_shelf[s * m_n_basins + b] = b;
502  }
503  }
504  }
505  }
506 
507  GlobalSum(m_grid->com, n_shelf_cells.data(),
508  n_shelf_cellsr.data(), m_n_shelves);
509  GlobalSum(m_grid->com, n_shelf_cells_per_basin.data(),
510  n_shelf_cells_per_basinr.data(), m_n_shelves*m_n_basins);
511  GlobalSum(m_grid->com, cfs_in_basins_per_shelf.data(),
512  cfs_in_basins_per_shelfr.data(), m_n_shelves*m_n_basins);
513  // copy data
514  n_shelf_cells = n_shelf_cellsr;
515  n_shelf_cells_per_basin = n_shelf_cells_per_basinr;
516  cfs_in_basins_per_shelf = cfs_in_basins_per_shelfr;
517 
518  for (int s = 0; s < m_n_shelves; s++) {
519  for (int b = 0; b < m_n_basins; b++) {
520  int sb = s * m_n_basins + b;
521  // remove ice shelf parts from the count that do not have a calving front in that basin
522  if (n_shelf_cells_per_basin[sb] > 0 and cfs_in_basins_per_shelf[sb] == 0) {
523  n_shelf_cells[s] -= n_shelf_cells_per_basin[sb];
524  n_shelf_cells_per_basin[sb] = 0;
525  }
526  }
527  }
528  }
529 
530  // now set potential temperature and salinity box 0:
531 
532  int low_temperature_counter = 0;
533  for (Points p(*m_grid); p; p.next()) {
534  const int i = p.i(), j = p.j();
535 
536  // make sure all temperatures are zero at the beginning of each time step
537  Toc_box0(i, j) = 0.0; // in K
538  Soc_box0(i, j) = 0.0; // in psu
539 
540  int s = shelf_mask.as_int(i, j);
541 
542  if (mask.as_int(i, j) == MASK_FLOATING and s > 0) {
543  // note: shelf_mask = 0 in lakes
544 
545  assert(n_shelf_cells[s] > 0);
546  double N = std::max(n_shelf_cells[s], 1); // protect from division by zero
547 
548  // weighted input depending on the number of shelf cells in each basin
549  for (int b = 1; b < m_n_basins; b++) { //Note: b=0 yields nan
550  int sb = s * m_n_basins + b;
551  Toc_box0(i, j) += basin_temperature[b] * n_shelf_cells_per_basin[sb] / N;
552  Soc_box0(i, j) += basin_salinity[b] * n_shelf_cells_per_basin[sb] / N;
553  }
554 
555  double theta_pm = physics.theta_pm(Soc_box0(i, j), physics.pressure(ice_thickness(i, j)));
556 
557  // temperature input for grounding line box should not be below pressure melting point
558  if (Toc_box0(i, j) < theta_pm) {
559  const double eps = 0.001;
560  // Setting Toc_box0 a little higher than theta_pm ensures that later equations are
561  // well solvable.
562  Toc_box0(i, j) = theta_pm + eps;
563  low_temperature_counter += 1;
564  }
565  }
566  }
567 
568  low_temperature_counter = GlobalSum(m_grid->com, low_temperature_counter);
569  if (low_temperature_counter > 0) {
570  m_log->message(2, "PICO ocean warning: temperature has been below pressure melting temperature in %d cases,\n"
571  "setting it to pressure melting temperature\n",
572  low_temperature_counter);
573  }
574 }
575 
576 /*!
577  * Use the simpler parameterization due to [@ref BeckmannGoosse2003] to set default
578  * sub-shelf temperature and melt rate values.
579  *
580  * At grid points containing floating ice not connected to the ocean, set the basal melt
581  * rate to zero and set basal temperature to the pressure melting point.
582  */
583 void Pico::beckmann_goosse(const PicoPhysics &physics,
584  const IceModelVec2S &ice_thickness,
585  const IceModelVec2Int &shelf_mask,
586  const IceModelVec2CellType &cell_type,
587  const IceModelVec2S &Toc_box0,
588  const IceModelVec2S &Soc_box0,
589  IceModelVec2S &basal_melt_rate,
590  IceModelVec2S &basal_temperature,
591  IceModelVec2S &Toc,
592  IceModelVec2S &Soc) {
593 
594  const double T0 = m_config->get_number("constants.fresh_water.melting_point_temperature"),
595  beta_CC = m_config->get_number("constants.ice.beta_Clausius_Clapeyron"),
596  g = m_config->get_number("constants.standard_gravity"),
597  ice_density = m_config->get_number("constants.ice.density");
598 
599  IceModelVec::AccessList list{ &ice_thickness, &cell_type, &shelf_mask, &Toc_box0, &Soc_box0,
600  &Toc, &Soc, &basal_melt_rate, &basal_temperature };
601 
602  for (Points p(*m_grid); p; p.next()) {
603  const int i = p.i(), j = p.j();
604 
605  if (cell_type.floating_ice(i, j)) {
606  if (shelf_mask.as_int(i, j) > 0) {
607  double pressure = physics.pressure(ice_thickness(i, j));
608 
609  basal_melt_rate(i, j) =
610  physics.melt_rate_beckmann_goosse(physics.theta_pm(Soc_box0(i, j), pressure), Toc_box0(i, j));
611  basal_temperature(i, j) = physics.T_pm(Soc_box0(i, j), pressure);
612 
613  // diagnostic outputs
614  Toc(i, j) = Toc_box0(i, j); // in Kelvin
615  Soc(i, j) = Soc_box0(i, j); // in psu
616  } else {
617  // Floating ice cells not connected to the ocean.
618  const double pressure = ice_density * g * ice_thickness(i, j); // FIXME issue #15
619 
620  basal_temperature(i, j) = T0 - beta_CC * pressure;
621  basal_melt_rate(i, j) = 0.0;
622  }
623  }
624  }
625 }
626 
627 
628 void Pico::process_box1(const PicoPhysics &physics,
629  const IceModelVec2S &ice_thickness,
630  const IceModelVec2Int &shelf_mask,
631  const IceModelVec2Int &box_mask,
632  const IceModelVec2S &Toc_box0,
633  const IceModelVec2S &Soc_box0,
634  IceModelVec2S &basal_melt_rate,
635  IceModelVec2S &basal_temperature,
636  IceModelVec2S &T_star,
637  IceModelVec2S &Toc,
638  IceModelVec2S &Soc,
639  IceModelVec2S &overturning) {
640 
641  std::vector<double> box1_area(m_n_shelves);
642 
643  compute_box_area(1, shelf_mask, box_mask, box1_area);
644 
645  IceModelVec::AccessList list{ &ice_thickness, &shelf_mask, &box_mask, &T_star, &Toc_box0, &Toc,
646  &Soc_box0, &Soc, &overturning, &basal_melt_rate, &basal_temperature };
647 
648  int n_Toc_failures = 0;
649 
650  // basal melt rate, ambient temperature and salinity and overturning calculation
651  // for each box1 grid cell.
652  for (Points p(*m_grid); p; p.next()) {
653  const int i = p.i(), j = p.j();
654 
655  int shelf_id = shelf_mask.as_int(i, j);
656 
657  if (box_mask.as_int(i, j) == 1 and shelf_id > 0) {
658 
659  const double pressure = physics.pressure(ice_thickness(i, j));
660 
661  T_star(i, j) = physics.T_star(Soc_box0(i, j), Toc_box0(i, j), pressure);
662 
663  auto Toc_box1 = physics.Toc_box1(box1_area[shelf_id], T_star(i, j), Soc_box0(i, j), Toc_box0(i, j));
664 
665  // This can only happen if T_star > 0.25*p_coeff, in particular T_star > 0
666  // which can only happen for values of Toc_box0 close to the local pressure melting point
667  if (Toc_box1.failed) {
668  m_log->message(5, "PICO ocean WARNING: negative square root argument at %d, %d\n"
669  "probably because of positive T_star=%f \n"
670  "Not aborting, but setting square root to 0... \n",
671  i, j, T_star(i, j));
672 
673  n_Toc_failures += 1;
674  }
675 
676  Toc(i, j) = Toc_box1.value;
677  Soc(i, j) = physics.Soc_box1(Toc_box0(i, j), Soc_box0(i, j), Toc(i, j)); // in psu
678 
679  overturning(i, j) = physics.overturning(Soc_box0(i, j), Soc(i, j), Toc_box0(i, j), Toc(i, j));
680 
681  // main outputs
682  basal_melt_rate(i, j) = physics.melt_rate(physics.theta_pm(Soc(i, j), pressure), Toc(i, j));
683  basal_temperature(i, j) = physics.T_pm(Soc(i, j), pressure);
684  }
685  }
686 
687  n_Toc_failures = GlobalSum(m_grid->com, n_Toc_failures);
688  if (n_Toc_failures > 0) {
689  m_log->message(2, "PICO ocean warning: square-root argument for temperature calculation "
690  "has been negative in %d cases!\n",
691  n_Toc_failures);
692  }
693 }
694 
696  const IceModelVec2S &ice_thickness,
697  const IceModelVec2Int &shelf_mask,
698  const IceModelVec2Int &box_mask,
699  IceModelVec2S &basal_melt_rate,
700  IceModelVec2S &basal_temperature,
701  IceModelVec2S &T_star,
702  IceModelVec2S &Toc,
703  IceModelVec2S &Soc) const {
704 
705  std::vector<double> overturning(m_n_shelves, 0.0);
706  std::vector<double> salinity(m_n_shelves, 0.0);
707  std::vector<double> temperature(m_n_shelves, 0.0);
708 
709  // get average overturning from box 1 that is used as input later
710  compute_box_average(1, m_overturning, shelf_mask, box_mask, overturning);
711 
712  std::vector<bool> use_beckmann_goosse(m_n_shelves);
713 
714  IceModelVec::AccessList list{ &ice_thickness, &shelf_mask, &box_mask, &T_star, &Toc,
715  &Soc, &basal_melt_rate, &basal_temperature };
716 
717  // Iterate over all boxes i for i > 1
718  for (int box = 2; box <= m_n_boxes; ++box) {
719 
720  compute_box_average(box - 1, Toc, shelf_mask, box_mask, temperature);
721  compute_box_average(box - 1, Soc, shelf_mask, box_mask, salinity);
722 
723  // find all the shelves where we should fall back to the Beckmann-Goosse
724  // parameterization
725  for (int s = 1; s < m_n_shelves; ++s) {
726  use_beckmann_goosse[s] = (salinity[s] == 0.0 or
727  temperature[s] == 0.0 or
728  overturning[s] == 0.0);
729  }
730 
731  std::vector<double> box_area;
732  compute_box_area(box, shelf_mask, box_mask, box_area);
733 
734  int n_beckmann_goosse_cells = 0;
735 
736  for (Points p(*m_grid); p; p.next()) {
737  const int i = p.i(), j = p.j();
738 
739  int shelf_id = shelf_mask.as_int(i, j);
740 
741  if (box_mask.as_int(i, j) == box and shelf_id > 0) {
742 
743  if (use_beckmann_goosse[shelf_id]) {
744  n_beckmann_goosse_cells += 1;
745  continue;
746  }
747 
748  // Get the input from previous box
749  const double
750  S_previous = salinity[shelf_id],
751  T_previous = temperature[shelf_id],
752  overturning_box1 = overturning[shelf_id];
753 
754  {
755  double pressure = physics.pressure(ice_thickness(i, j));
756 
757  // diagnostic outputs
758  T_star(i, j) = physics.T_star(S_previous, T_previous, pressure);
759  Toc(i, j) = physics.Toc(box_area[shelf_id], T_previous, T_star(i, j), overturning_box1, S_previous);
760  Soc(i, j) = physics.Soc(S_previous, T_previous, Toc(i, j));
761 
762  // main outputs: basal melt rate and temperature
763  basal_melt_rate(i, j) = physics.melt_rate(physics.theta_pm(Soc(i, j), pressure), Toc(i, j));
764  basal_temperature(i, j) = physics.T_pm(Soc(i, j), pressure);
765  }
766  }
767  } // loop over grid points
768 
769  n_beckmann_goosse_cells = GlobalSum(m_grid->com, n_beckmann_goosse_cells);
770  if (n_beckmann_goosse_cells > 0) {
771  m_log->message(2, "PICO ocean WARNING: box %d, no boundary data from previous box in %d case(s)!\n"
772  "switching to Beckmann Goosse (2003) meltrate calculation\n",
773  box, n_beckmann_goosse_cells);
774  }
775 
776  } // loop over boxes
777 }
778 
779 
780 // Write diagnostic variables to extra files if requested
782 
783  DiagnosticList result = {
784  { "basins", Diagnostic::wrap(m_geometry.basin_mask()) },
785  { "pico_overturning", Diagnostic::wrap(m_overturning) },
786  { "pico_salinity_box0", Diagnostic::wrap(m_Soc_box0) },
787  { "pico_temperature_box0", Diagnostic::wrap(m_Toc_box0) },
788  { "pico_box_mask", Diagnostic::wrap(m_geometry.box_mask()) },
789  { "pico_shelf_mask", Diagnostic::wrap(m_geometry.ice_shelf_mask()) },
790  { "pico_ice_rise_mask", Diagnostic::wrap(m_geometry.ice_rise_mask()) },
791  { "pico_basal_melt_rate", Diagnostic::wrap(m_basal_melt_rate) },
792  { "pico_contshelf_mask", Diagnostic::wrap(m_geometry.continental_shelf_mask()) },
793  { "pico_salinity", Diagnostic::wrap(m_Soc) },
794  { "pico_temperature", Diagnostic::wrap(m_Toc) },
795  { "pico_T_star", Diagnostic::wrap(m_T_star) },
796  { "pico_basal_temperature", Diagnostic::wrap(*m_shelf_base_temperature) },
797  };
798 
799  return combine(result, OceanModel::diagnostics_impl());
800 }
801 
802 /*!
803  * For each shelf, compute average of a given field over the box with id `box_id`.
804  *
805  * This method is used to get inputs from a previous box for the next one.
806  */
808  const IceModelVec2S &field,
809  const IceModelVec2Int &shelf_mask,
810  const IceModelVec2Int &box_mask,
811  std::vector<double> &result) const {
812 
813  IceModelVec::AccessList list{ &field, &shelf_mask, &box_mask };
814 
815  // fill results with zeros
816  result.resize(m_n_shelves);
817  for (int s = 0; s < m_n_shelves; ++s) {
818  result[s] = 0.0;
819  }
820 
821  // compute the sum of field in each shelf's box box_id
822  std::vector<int> n_cells(m_n_shelves);
823  {
824  std::vector<int> n_cells_per_box(m_n_shelves, 0);
825  for (Points p(*m_grid); p; p.next()) {
826  const int i = p.i(), j = p.j();
827 
828  int shelf_id = shelf_mask.as_int(i, j);
829 
830  if (box_mask.as_int(i, j) == box_id) {
831  n_cells_per_box[shelf_id] += 1;
832  result[shelf_id] += field(i, j);
833  }
834  }
835  GlobalSum(m_grid->com, n_cells_per_box.data(), n_cells.data(), m_n_shelves);
836  }
837 
838  {
839  std::vector<double> tmp(m_n_shelves);
840  GlobalSum(m_grid->com, result.data(), tmp.data(), m_n_shelves);
841  // copy data
842  result = tmp;
843  }
844 
845  for (int s = 0; s < m_n_shelves; ++s) {
846  if (n_cells[s] > 0) {
847  result[s] /= static_cast<double>(n_cells[s]);
848  }
849  }
850 }
851 
852 /*!
853  * For all shelves compute areas of boxes with id `box_id`.
854  *
855  * @param[in] box_mask box index mask
856  * @param[out] result resulting box areas.
857  *
858  * Note: shelf and box indexes start from 1.
859  */
860 void Pico::compute_box_area(int box_id,
861  const IceModelVec2Int &shelf_mask,
862  const IceModelVec2Int &box_mask,
863  std::vector<double> &result) const {
864  result.resize(m_n_shelves);
865  IceModelVec::AccessList list{ &shelf_mask, &box_mask };
866 
867  auto cell_area = m_grid->cell_area();
868 
869  for (Points p(*m_grid); p; p.next()) {
870  const int i = p.i(), j = p.j();
871 
872  int shelf_id = shelf_mask.as_int(i, j);
873 
874  if (shelf_id > 0 and box_mask.as_int(i, j) == box_id) {
875  result[shelf_id] += cell_area;
876  }
877  }
878 
879  // compute GlobalSum from index 1 to index m_n_shelves-1
880  std::vector<double> result1(m_n_shelves);
881  GlobalSum(m_grid->com, &result[1], &result1[1], m_n_shelves-1);
882  // copy data
883  for (int i = 1; i < m_n_shelves; i++) {
884  result[i] = result1[i];
885  }
886 }
887 
888 } // end of namespace ocean
889 } // end of namespace pism
Makes sure that we call begin_access() and end_access() for all accessed IceModelVecs.
Definition: iceModelVec.hh:59
const Config::ConstPtr m_config
configuration database used by this component
Definition: Component.hh:138
const Logger::ConstPtr m_log
logger (for easy access)
Definition: Component.hh:142
const IceGrid::ConstPtr m_grid
grid used by this component
Definition: Component.hh:136
static Ptr wrap(const T &input)
Definition: Diagnostic.hh:155
High-level PISM I/O class.
Definition: File.hh:51
IceModelVec2CellType cell_type
Definition: Geometry.hh:57
IceModelVec2S bed_elevation
Definition: Geometry.hh:49
IceModelVec2S ice_thickness
Definition: Geometry.hh:53
std::shared_ptr< const IceGrid > ConstPtr
Definition: IceGrid.hh:233
bool floating_ice(int i, int j) const
"Cell type" mask. Adds convenience methods to IceModelVec2Int.
stencils::Box< int > box(int i, int j) const
int as_int(int i, int j) const
stencils::Star< int > star(int i, int j) const
A simple class "hiding" the fact that the mask is stored as floating-point scalars (instead of intege...
Definition: iceModelVec.hh:389
stencils::Box< double > box(int i, int j) const
static std::shared_ptr< IceModelVec2T > ForcingField(IceGrid::ConstPtr grid, const File &file, const std::string &short_name, const std::string &standard_name, int max_buffer_size, bool periodic, InterpolationType interpolation_type=PIECEWISE_CONSTANT)
void update_ghosts()
Updates ghost points.
Definition: iceModelVec.cc:669
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition: iceModelVec.cc:533
void set_attrs(const std::string &pism_intent, const std::string &long_name, const std::string &units, const std::string &glaciological_units, const std::string &standard_name, unsigned int component)
Sets NetCDF attributes of an IceModelVec object.
Definition: iceModelVec.cc:399
void set(double c)
Result: v[j] <- c for all j.
Definition: iceModelVec.cc:683
IceGrid::ConstPtr grid() const
Definition: iceModelVec.cc:128
void define(const File &file, IO_Type default_type=PISM_DOUBLE) const
Define variables corresponding to an IceModelVec in a file opened using file.
Definition: iceModelVec.cc:523
void write(const std::string &filename) const
Definition: iceModelVec.cc:822
Combines the max. time step with the flag indicating if a restriction is active. Makes is possible to...
Definition: MaxTimestep.hh:31
double get_number(const std::string &name) const
Get a single-valued scalar attribute.
IceModelVec2S::Ptr m_shelf_base_mass_flux
IceModelVec2S::Ptr m_shelf_base_temperature
virtual void define_model_state_impl(const File &output) const
The default (empty implementation).
Definition: OceanModel.cc:128
IceModelVec2S::Ptr m_water_column_pressure
Definition: OceanModel.hh:70
virtual DiagnosticList diagnostics_impl() const
Definition: OceanModel.cc:224
static void compute_average_water_column_pressure(const Geometry &geometry, double ice_density, double water_density, double g, IceModelVec2S &result)
Definition: OceanModel.cc:246
virtual void write_model_state_impl(const File &output) const
The default (empty implementation).
Definition: OceanModel.cc:136
A very rudimentary PISM ocean model.
Definition: OceanModel.hh:36
const IceModelVec2Int & continental_shelf_mask() const
Definition: PicoGeometry.cc:58
const IceModelVec2Int & basin_mask() const
Definition: PicoGeometry.cc:74
const IceModelVec2Int & ice_shelf_mask() const
Definition: PicoGeometry.cc:66
void update(const IceModelVec2S &bed_elevation, const IceModelVec2CellType &cell_type)
Definition: PicoGeometry.cc:93
const IceModelVec2Int & box_mask() const
Definition: PicoGeometry.cc:62
const IceModelVec2Int & ice_rise_mask() const
Definition: PicoGeometry.cc:70
double pressure(double ice_thickness) const
Definition: PicoPhysics.cc:90
double overturning_coeff() const
Definition: PicoPhysics.cc:195
double ice_density() const
Definition: PicoPhysics.cc:207
double S_dummy() const
Definition: PicoPhysics.cc:203
TocBox1 Toc_box1(double area, double T_star, double Soc_box0, double Toc_box0) const
Definition: PicoPhysics.cc:95
double Toc(double box_area, double temperature, double T_star, double overturning, double salinity) const
Definition: PicoPhysics.cc:120
double gamma_T() const
Definition: PicoPhysics.cc:191
double continental_shelf_depth() const
Definition: PicoPhysics.cc:211
double melt_rate(double pm_point, double Toc) const
equation 8 in the PICO paper.
Definition: PicoPhysics.cc:155
double overturning(double Soc_box0, double Soc, double Toc_box0, double Toc) const
equation 3 in the PICO paper. See also equation 4.
Definition: PicoPhysics.cc:169
double theta_pm(double salinity, double pressure) const
Definition: PicoPhysics.cc:142
double Soc_box1(double Toc_box0, double Soc_box0, double Toc) const
Definition: PicoPhysics.cc:129
double T_pm(double salinity, double pressure) const
Definition: PicoPhysics.cc:149
double melt_rate_beckmann_goosse(double pot_pm_point, double Toc) const
Beckmann & Goosse meltrate.
Definition: PicoPhysics.cc:161
double T_star(double salinity, double temperature, double pressure) const
See equation A6 and lines before in PICO paper.
Definition: PicoPhysics.cc:175
double T_dummy() const
Definition: PicoPhysics.cc:199
double Soc(double salinity, double temperature, double Toc) const
Definition: PicoPhysics.cc:134
std::shared_ptr< IceModelVec2T > m_salinity_ocean
Definition: Pico.hh:64
void process_box1(const PicoPhysics &physics, const IceModelVec2S &ice_thickness, const IceModelVec2Int &shelf_mask, const IceModelVec2Int &box_mask, const IceModelVec2S &Toc_box0, const IceModelVec2S &Soc_box0, IceModelVec2S &basal_melt_rate, IceModelVec2S &basal_temperature, IceModelVec2S &T_star, IceModelVec2S &Toc, IceModelVec2S &Soc, IceModelVec2S &overturning)
Definition: Pico.cc:628
IceModelVec2S m_Toc_box0
Definition: Pico.hh:58
void compute_box_area(int box_id, const IceModelVec2Int &shelf_mask, const IceModelVec2Int &box_mask, std::vector< double > &result) const
Definition: Pico.cc:860
Pico(IceGrid::ConstPtr g)
Definition: Pico.cc:51
void write_model_state_impl(const File &output) const
The default (empty implementation).
Definition: Pico.cc:189
MaxTimestep max_timestep_impl(double t) const
Definition: Pico.cc:366
void set_ocean_input_fields(const PicoPhysics &physics, const IceModelVec2S &ice_thickness, const IceModelVec2CellType &mask, const IceModelVec2Int &basin_mask, const IceModelVec2Int &shelf_mask, const std::vector< double > &basin_temperature, const std::vector< double > &basin_salinity, IceModelVec2S &Toc_box0, IceModelVec2S &Soc_box0) const
Set ocean ocean input from box 0 as boundary condition for box 1.
Definition: Pico.cc:463
void beckmann_goosse(const PicoPhysics &physics, const IceModelVec2S &ice_thickness, const IceModelVec2Int &shelf_mask, const IceModelVec2CellType &cell_type, const IceModelVec2S &Toc_box0, const IceModelVec2S &Soc_box0, IceModelVec2S &basal_melt_rate, IceModelVec2S &basal_temperature, IceModelVec2S &Toc, IceModelVec2S &Soc)
Definition: Pico.cc:583
IceModelVec2S m_overturning
Definition: Pico.hh:59
void define_model_state_impl(const File &output) const
The default (empty implementation).
Definition: Pico.cc:179
void compute_box_average(int box_id, const IceModelVec2S &field, const IceModelVec2Int &shelf_mask, const IceModelVec2Int &box_mask, std::vector< double > &result) const
Definition: Pico.cc:807
void update_impl(const Geometry &geometry, double t, double dt)
Definition: Pico.cc:250
std::shared_ptr< IceModelVec2T > m_theta_ocean
Definition: Pico.hh:64
IceModelVec2S m_basal_melt_rate
Definition: Pico.hh:60
PicoGeometry m_geometry
Definition: Pico.hh:62
void init_impl(const Geometry &geometry)
Definition: Pico.cc:134
IceModelVec2S m_Toc
Definition: Pico.hh:58
std::map< std::string, Diagnostic::Ptr > diagnostics_impl() const
Definition: Pico.cc:781
IceModelVec2S m_Soc_box0
Definition: Pico.hh:57
IceModelVec2S m_T_star
Definition: Pico.hh:58
void process_other_boxes(const PicoPhysics &physics, const IceModelVec2S &ice_thickness, const IceModelVec2Int &shelf_mask, const IceModelVec2Int &box_mask, IceModelVec2S &basal_melt_rate, IceModelVec2S &basal_temperature, IceModelVec2S &T_star, IceModelVec2S &Toc, IceModelVec2S &Soc) const
Definition: Pico.cc:695
IceModelVec2S m_Soc
Definition: Pico.hh:57
void compute_ocean_input_per_basin(const PicoPhysics &physics, const IceModelVec2Int &basin_mask, const IceModelVec2Int &continental_shelf_mask, const IceModelVec2S &salinity_ocean, const IceModelVec2S &theta_ocean, std::vector< double > &temperature, std::vector< double > &salinity) const
Compute temperature and salinity input from ocean data by averaging.
Definition: Pico.cc:379
static const double beta_CC
Definition: exactTestO.c:23
bool ocean(int M)
An ocean cell (floating ice or ice-free).
Definition: Mask.hh:39
static void extend_basal_melt_rates(const IceModelVec2CellType &cell_type, IceModelVec2S &basal_melt_rate)
Definition: Pico.cc:204
double max(const IceModelVec2S &input)
Finds maximum over all the values in an IceModelVec2S object. Ignores ghosts.
static const double g
Definition: exactTestP.cc:39
std::map< std::string, Diagnostic::Ptr > DiagnosticList
Definition: Diagnostic.hh:117
@ MASK_FLOATING
Definition: Mask.hh:33
@ MASK_ICE_FREE_OCEAN
Definition: Mask.hh:34
@ MASK_GROUNDED
Definition: Mask.hh:32
@ PISM_NETCDF3
Definition: IO_Flags.hh:41
double min(const IceModelVec2S &input)
Finds minimum over all the values in an IceModelVec2S object. Ignores ghosts.
@ PISM_READONLY
open an existing file for reading only
Definition: IO_Flags.hh:49
T combine(const T &a, const T &b)
void GlobalSum(MPI_Comm comm, double *local, double *result, int count)
@ WITHOUT_GHOSTS
Definition: iceModelVec.hh:49
@ WITH_GHOSTS
Definition: iceModelVec.hh:49
std::string filename
Definition: options.hh:33
int count
Definition: test_cube.c:16