33 #include <gsl/gsl_math.h>
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"
59 m_basal_melt_rate(m_grid,
"pico_basal_melt_rate",
WITH_GHOSTS),
68 auto buffer_size =
static_cast<int>(
m_config->get_number(
"input.forcing.buffer_size"));
90 "potential temperature of the adjacent ocean",
91 "Kelvin",
"Kelvin",
"", 0);
94 "salinity of the adjacent ocean",
95 "g/kg",
"g/kg",
"", 0);
99 "g/kg",
"g/kg",
"ocean salinity field", 0);
104 "g/kg",
"g/kg",
"", 0);
118 "degree C",
"degree C",
"", 0);
122 "m^3 s-1",
"m^3 s-1",
"", 0);
126 "m s-1",
"m year-1",
"", 0);
131 m_n_boxes =
static_cast<int>(
m_config->get_number(
"ocean.pico.number_of_boxes"));
137 m_log->message(2,
"* Initializing the Potsdam Ice-shelf Cavity mOdel for the ocean ...\n");
150 m_log->message(4,
"PICO basin min=%f, max=%f\n",
156 m_log->message(2,
" -Using %d drainage basins and values: \n"
157 " gamma_T= %.2e, overturning_coeff = %.2e... \n",
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",
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");
207 auto grid = basal_melt_rate.
grid();
217 const int i = p.i(), j = p.j();
219 auto M = cell_type.
box(i, j);
221 bool potential_partially_filled_cell =
228 if (potential_partially_filled_cell) {
229 auto BMR = basal_melt_rate.
box(i, j);
232 double melt_sum = 0.0;
244 basal_melt_rate(i, j) = melt_sum / N;
260 double T_fill_value =
m_config->get_number(
"constants.fresh_water.melting_point_temperature");
291 std::vector<double> basin_temperature(
m_n_basins);
292 std::vector<double> basin_salinity(
m_n_basins);
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");
384 std::vector<double> &temperature,
385 std::vector<double> &salinity)
const {
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;
405 const int i = p.i(), j = p.j();
407 if (continental_shelf_mask.
as_int(i, j) == 2) {
408 int basin_id = basin_mask.
as_int(i, j);
410 count[basin_id] += 1;
411 salinity[basin_id] += salinity_ocean(i, j);
412 temperature[basin_id] += theta_ocean(i, j);
426 salinity = salinityr;
427 temperature = temperaturer;
431 temperature[0] = physics.
T_dummy();
432 salinity[0] = physics.
S_dummy();
435 for (
int basin_id = 1; basin_id <
m_n_basins; basin_id++) {
437 if (
count[basin_id] != 0) {
438 salinity[basin_id] /=
count[basin_id];
439 temperature[basin_id] /=
count[basin_id];
441 m_log->message(5,
" %d: temp =%.3f, salinity=%.3f\n", basin_id, temperature[basin_id], salinity[basin_id]);
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",
450 temperature[basin_id] = physics.
T_dummy();
451 salinity[basin_id] = physics.
S_dummy();
468 const std::vector<double> &basin_temperature,
469 const std::vector<double> &basin_salinity,
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);
495 auto M = mask.
star(i, j);
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;
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;
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;
532 int low_temperature_counter = 0;
534 const int i = p.i(), j = p.j();
537 Toc_box0(i, j) = 0.0;
538 Soc_box0(i, j) = 0.0;
540 int s = shelf_mask.
as_int(i, j);
545 assert(n_shelf_cells[s] > 0);
546 double N =
std::max(n_shelf_cells[s], 1);
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;
555 double theta_pm = physics.
theta_pm(Soc_box0(i, j), physics.
pressure(ice_thickness(i, j)));
558 if (Toc_box0(i, j) < theta_pm) {
559 const double eps = 0.001;
562 Toc_box0(i, j) = theta_pm + eps;
563 low_temperature_counter += 1;
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);
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");
600 &Toc, &Soc, &basal_melt_rate, &basal_temperature };
603 const int i = p.i(), j = p.j();
606 if (shelf_mask.
as_int(i, j) > 0) {
607 double pressure = physics.
pressure(ice_thickness(i, j));
609 basal_melt_rate(i, j) =
611 basal_temperature(i, j) = physics.
T_pm(Soc_box0(i, j), pressure);
614 Toc(i, j) = Toc_box0(i, j);
615 Soc(i, j) = Soc_box0(i, j);
618 const double pressure = ice_density *
g * ice_thickness(i, j);
620 basal_temperature(i, j) = T0 -
beta_CC * pressure;
621 basal_melt_rate(i, j) = 0.0;
646 &Soc_box0, &Soc, &overturning, &basal_melt_rate, &basal_temperature };
648 int n_Toc_failures = 0;
653 const int i = p.i(), j = p.j();
655 int shelf_id = shelf_mask.
as_int(i, j);
657 if (box_mask.
as_int(i, j) == 1 and shelf_id > 0) {
659 const double pressure = physics.
pressure(ice_thickness(i, j));
661 T_star(i, j) = physics.
T_star(Soc_box0(i, j), Toc_box0(i, j), pressure);
663 auto Toc_box1 = physics.
Toc_box1(box1_area[shelf_id], T_star(i, j), Soc_box0(i, j), Toc_box0(i, j));
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",
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));
679 overturning(i, j) = physics.
overturning(Soc_box0(i, j), Soc(i, j), Toc_box0(i, j), Toc(i, j));
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);
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",
712 std::vector<bool> use_beckmann_goosse(
m_n_shelves);
715 &Soc, &basal_melt_rate, &basal_temperature };
718 for (
int box = 2; box <=
m_n_boxes; ++box) {
726 use_beckmann_goosse[s] = (salinity[s] == 0.0 or
727 temperature[s] == 0.0 or
728 overturning[s] == 0.0);
731 std::vector<double> box_area;
734 int n_beckmann_goosse_cells = 0;
737 const int i = p.i(), j = p.j();
739 int shelf_id = shelf_mask.
as_int(i, j);
741 if (box_mask.
as_int(i, j) == box and shelf_id > 0) {
743 if (use_beckmann_goosse[shelf_id]) {
744 n_beckmann_goosse_cells += 1;
750 S_previous = salinity[shelf_id],
751 T_previous = temperature[shelf_id],
752 overturning_box1 = overturning[shelf_id];
755 double pressure = physics.
pressure(ice_thickness(i, j));
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));
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);
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);
811 std::vector<double> &result)
const {
826 const int i = p.i(), j = p.j();
828 int shelf_id = shelf_mask.
as_int(i, j);
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);
846 if (n_cells[s] > 0) {
847 result[s] /=
static_cast<double>(n_cells[s]);
863 std::vector<double> &result)
const {
867 auto cell_area =
m_grid->cell_area();
870 const int i = p.i(), j = p.j();
872 int shelf_id = shelf_mask.
as_int(i, j);
874 if (shelf_id > 0 and box_mask.
as_int(i, j) == box_id) {
875 result[shelf_id] += cell_area;
884 result[i] = result1[i];
Makes sure that we call begin_access() and end_access() for all accessed IceModelVecs.
const Config::ConstPtr m_config
configuration database used by this component
const Logger::ConstPtr m_log
logger (for easy access)
const IceGrid::ConstPtr m_grid
grid used by this component
static Ptr wrap(const T &input)
High-level PISM I/O class.
IceModelVec2CellType cell_type
IceModelVec2S bed_elevation
IceModelVec2S ice_thickness
std::shared_ptr< const IceGrid > ConstPtr
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...
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.
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
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.
void set(double c)
Result: v[j] <- c for all j.
IceGrid::ConstPtr grid() const
void define(const File &file, IO_Type default_type=PISM_DOUBLE) const
Define variables corresponding to an IceModelVec in a file opened using file.
void write(const std::string &filename) const
Combines the max. time step with the flag indicating if a restriction is active. Makes is possible to...
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).
IceModelVec2S::Ptr m_water_column_pressure
virtual DiagnosticList diagnostics_impl() const
static void compute_average_water_column_pressure(const Geometry &geometry, double ice_density, double water_density, double g, IceModelVec2S &result)
virtual void write_model_state_impl(const File &output) const
The default (empty implementation).
A very rudimentary PISM ocean model.
const IceModelVec2Int & continental_shelf_mask() const
const IceModelVec2Int & basin_mask() const
const IceModelVec2Int & ice_shelf_mask() const
void update(const IceModelVec2S &bed_elevation, const IceModelVec2CellType &cell_type)
const IceModelVec2Int & box_mask() const
const IceModelVec2Int & ice_rise_mask() const
double pressure(double ice_thickness) const
double overturning_coeff() const
double ice_density() const
TocBox1 Toc_box1(double area, double T_star, double Soc_box0, double Toc_box0) const
double Toc(double box_area, double temperature, double T_star, double overturning, double salinity) const
double continental_shelf_depth() const
double melt_rate(double pm_point, double Toc) const
equation 8 in the PICO paper.
double overturning(double Soc_box0, double Soc, double Toc_box0, double Toc) const
equation 3 in the PICO paper. See also equation 4.
double theta_pm(double salinity, double pressure) const
double Soc_box1(double Toc_box0, double Soc_box0, double Toc) const
double T_pm(double salinity, double pressure) const
double melt_rate_beckmann_goosse(double pot_pm_point, double Toc) const
Beckmann & Goosse meltrate.
double T_star(double salinity, double temperature, double pressure) const
See equation A6 and lines before in PICO paper.
double Soc(double salinity, double temperature, double Toc) const
std::shared_ptr< IceModelVec2T > m_salinity_ocean
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)
void compute_box_area(int box_id, const IceModelVec2Int &shelf_mask, const IceModelVec2Int &box_mask, std::vector< double > &result) const
Pico(IceGrid::ConstPtr g)
void write_model_state_impl(const File &output) const
The default (empty implementation).
MaxTimestep max_timestep_impl(double t) const
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.
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)
IceModelVec2S m_overturning
void define_model_state_impl(const File &output) const
The default (empty implementation).
void compute_box_average(int box_id, const IceModelVec2S &field, const IceModelVec2Int &shelf_mask, const IceModelVec2Int &box_mask, std::vector< double > &result) const
void update_impl(const Geometry &geometry, double t, double dt)
std::shared_ptr< IceModelVec2T > m_theta_ocean
IceModelVec2S m_basal_melt_rate
void init_impl(const Geometry &geometry)
std::map< std::string, Diagnostic::Ptr > diagnostics_impl() const
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
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.
static const double beta_CC
bool ocean(int M)
An ocean cell (floating ice or ice-free).
static void extend_basal_melt_rates(const IceModelVec2CellType &cell_type, IceModelVec2S &basal_melt_rate)
double max(const IceModelVec2S &input)
Finds maximum over all the values in an IceModelVec2S object. Ignores ghosts.
std::map< std::string, Diagnostic::Ptr > DiagnosticList
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
T combine(const T &a, const T &b)
void GlobalSum(MPI_Comm comm, double *local, double *result, int count)