23 #include "pism/util/IceGrid.hh"
24 #include "pism/util/pism_options.hh"
25 #include "pism/util/Vars.hh"
26 #include "pism/util/Time.hh"
27 #include "pism/coupler/AtmosphereModel.hh"
28 #include "pism/util/Mask.hh"
29 #include "pism/util/io/File.hh"
31 #include "pism/util/error_handling.hh"
32 #include "pism/util/io/io_helpers.hh"
33 #include "pism/util/pism_utilities.hh"
34 #include "pism/util/IceModelVec2CellType.hh"
35 #include "pism/geometry/Geometry.hh"
36 #include "pism/util/Context.hh"
44 std::shared_ptr<atmosphere::AtmosphereModel> input)
58 bool use_fausto_params =
m_config->get_flag(
"surface.pdd.fausto.enabled");
60 std::string method =
m_config->get_string(
"surface.pdd.method");
62 if (method ==
"repeatable_random_process") {
64 }
else if (method ==
"random_process") {
70 if (use_fausto_params) {
75 std::string sd_file =
m_config->get_string(
"surface.pdd.std_dev.file");
77 if (not sd_file.empty()) {
78 bool sd_periodic =
m_config->get_flag(
"surface.pdd.std_dev.periodic");
79 int max_buffer_size = (
unsigned int)
m_config->get_number(
"input.forcing.buffer_size");
94 "standard deviation of near-surface air temperature",
95 "Kelvin",
"Kelvin",
"", 0);
98 "instantaneous surface mass balance (accumulation/ablation) rate",
99 "kg m-2 s-1",
"kg m-2 s-1",
100 "land_ice_surface_specific_mass_balance_flux", 0);
105 "snow cover depth (set to zero once a year)",
130 "* Initializing the default temperature-index, PDD-based surface processes scheme.\n"
131 " Precipitation and 2m air temperature provided by atmosphere are inputs.\n"
132 " Surface mass balance and ice upper surface temperature are outputs.\n"
133 " See PISM User's Manual for control of degree-day factors.\n");
136 " Computing number of positive degree-days by: %s.\n",
141 " Setting PDD parameters from [Faustoetal2009].\n");
144 " Using default PDD parameters.\n");
150 std::string sd_file =
m_config->get_string(
"surface.pdd.std_dev.file");
151 if (sd_file.empty()) {
153 " Using constant standard deviation of near-surface temperature.\n");
162 " Reading standard deviation of near-surface temperature from '%s'...\n",
165 bool sd_periodic =
m_config->get_flag(
"surface.pdd.std_dev.periodic");
173 std::string firn_file =
m_config->get_string(
"surface.pdd.firn_depth_file");
176 if (not firn_file.empty()) {
178 "surface.pdd.firn_depth_file is not allowed when"
179 " re-starting from a PISM output file.");
188 if (firn_file.empty()) {
197 if (firn_file.empty()) {
226 balance_year_start_day =
m_config->get_number(
"surface.pdd.balance_year_start_day"),
228 year_start =
m_grid->ctx()->time()->calendar_year_start(time),
229 balance_year_start = year_start + (balance_year_start_day - 1.0) * one_day;
231 if (balance_year_start > time) {
232 return balance_year_start;
234 return m_grid->ctx()->time()->increment_date(balance_year_start, 1);
249 int N =
m_mbscheme->get_timeseries_length(dt);
251 const double dtseries = dt / N;
252 std::vector<double> ts(N), T(N),
S(N), P(N), PDDs(N);
253 for (
int k = 0;
k < N; ++
k) {
254 ts[
k] = t +
k * dtseries;
271 sigmalapserate =
m_config->get_number(
"surface.pdd.std_dev.lapse_lat_rate"),
272 sigmabaselat =
m_config->get_number(
"surface.pdd.std_dev.lapse_lat_base");
275 if (fausto_greve or sigmalapserate != 0.0) {
286 fausto_greve->
update_temp_mj(*surface_altitude, *latitude, *longitude);
295 const double ice_density =
m_config->get_number(
"constants.ice.density");
300 const int i = p.i(), j = p.j();
307 for (
int k = 0;
k < N; ++
k) {
317 for (
int k = 0;
k < N; ++
k) {
318 P[
k] = P[
k] / ice_density;
326 double tmp = (*m_air_temp_sd)(i, j);
327 for (
int k = 0;
k < N; ++
k) {
339 if (sigmalapserate != 0.0) {
340 double lat = (*latitude)(i, j);
341 for (
int k = 0;
k < N; ++
k) {
342 S[
k] += sigmalapserate * (lat - sigmabaselat);
344 (*m_air_temp_sd)(i, j) =
S[0];
349 for (
int k = 0;
k < N; ++
k) {
355 (*m_air_temp_sd)(i, j) =
S[0];
362 for (
int k = 0;
k < N; ++
k) {
393 for (
int k = 0;
k < N; ++
k) {
394 if (ts[
k] >= next_snow_depth_reset) {
396 while (next_snow_depth_reset <= ts[
k]) {
397 next_snow_depth_reset =
m_grid->ctx()->time()->increment_date(next_snow_depth_reset, 1);
435 (*m_accumulation)(i, j) = A * ice_density;
436 (*m_melt)(i, j) = M * ice_density;
437 (*m_runoff)(i, j) = R * ice_density;
503 namespace diagnostics {
514 ?
"surface_melt_flux"
515 :
"surface_melt_rate",
522 name =
"surface_melt_flux",
523 long_name =
"surface melt, averaged over the reporting interval",
524 standard_name =
"surface_snow_and_ice_melt_flux",
525 accumulator_units =
"kg m-2",
526 internal_units =
"kg m-2 second-1",
527 external_units =
"kg m-2 year-1";
529 name =
"surface_melt_rate";
531 accumulator_units =
"kg",
532 internal_units =
"kg second-1";
533 external_units =
"Gt year-1" ;
539 set_attrs(long_name, standard_name, internal_units, external_units, 0);
540 m_vars[0][
"cell_methods"] =
"time: mean";
550 double cell_area =
m_grid->cell_area();
555 const int i = p.i(), j = p.j();
575 ?
"surface_runoff_flux"
576 :
"surface_runoff_rate",
582 name =
"surface_runoff_flux",
583 long_name =
"surface runoff, averaged over the reporting interval",
584 standard_name =
"surface_runoff_flux",
585 accumulator_units =
"kg m-2",
586 internal_units =
"kg m-2 second-1",
587 external_units =
"kg m-2 year-1";
589 name =
"surface_runoff_rate";
591 accumulator_units =
"kg",
592 internal_units =
"kg second-1";
593 external_units =
"Gt year-1" ;
599 set_attrs(long_name, standard_name, internal_units, external_units, 0);
600 m_vars[0][
"cell_methods"] =
"time: mean";
610 double cell_area =
m_grid->cell_area();
615 const int i = p.i(), j = p.j();
620 return runoff_amount;
635 ?
"surface_accumulation_flux"
636 :
"surface_accumulation_rate",
643 name =
"surface_accumulation_flux",
644 long_name =
"accumulation (precipitation minus rain), averaged over the reporting interval",
645 accumulator_units =
"kg m-2",
646 internal_units =
"kg m-2 second-1",
647 external_units =
"kg m-2 year-1";
649 name =
"surface_accumulation_rate";
650 accumulator_units =
"kg",
651 internal_units =
"kg second-1";
652 external_units =
"Gt year-1" ;
659 set_attrs(long_name,
"", internal_units, external_units, 0);
660 m_vars[0][
"cell_methods"] =
"time: mean";
670 double cell_area =
m_grid->cell_area();
675 const int i = p.i(), j = p.j();
680 return accumulation_amount;
696 double cell_area = grid->cell_area();
703 const int i = p.i(), j = p.j();
705 result += input(i, j) * cell_area;
720 m_variable[
"long_name"] =
"surface accumulation rate (PDD model)";
737 m_variable[
"long_name"] =
"surface melt rate (PDD model)";
754 m_variable[
"long_name"] =
"surface runoff rate (PDD model)";
788 {
"surface_accumulation_rate",
TSDiagnostic::Ptr(
new TotalSurfaceAccumulation(
this))},
Makes sure that we call begin_access() and end_access() for all accessed IceModelVecs.
const units::System::Ptr m_sys
unit system used by this component
const Config::ConstPtr m_config
configuration database used by this component
const Logger::ConstPtr m_log
logger (for easy access)
virtual void regrid(const std::string &module_name, IceModelVec &variable, RegriddingFlag flag=NO_REGRID_WITHOUT_REGRID_VARS)
const IceGrid::ConstPtr m_grid
grid used by this component
DiagnosticList diagnostics() const
IceModelVec2S m_accumulator
static Ptr wrap(const T &input)
double m_fill_value
fill value (used often enough to justify storing it)
const units::System::Ptr m_sys
the unit system
double to_internal(double x) const
std::vector< SpatialVariableMetadata > m_vars
metadata corresponding to NetCDF variables
std::shared_ptr< Diagnostic > Ptr
IceGrid::ConstPtr m_grid
the grid
void set_attrs(const std::string &long_name, const std::string &standard_name, const std::string &units, const std::string &glaciological_units, unsigned int N=0)
A method for setting common variable attributes.
High-level PISM I/O class.
IceModelVec2S ice_surface_elevation
IceModelVec2CellType cell_type
IceModelVec2S ice_thickness
std::shared_ptr< const IceGrid > ConstPtr
bool icy(int i, int j) const
bool ice_free_ocean(int i, int j) const
"Cell type" mask. Adds convenience methods to IceModelVec2Int.
void add(double alpha, const IceModelVec2S &x)
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)
static std::shared_ptr< IceModelVec2T > Constant(IceGrid::ConstPtr grid, const std::string &short_name, double value)
A class for storing and accessing 2D time-series (for climate forcing)
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 regrid(const std::string &filename, RegriddingFlag flag, double default_value=0.0)
void set(double c)
Result: v[j] <- c for all j.
void read(const std::string &filename, unsigned int time)
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...
void failed()
Indicates a failure of a parallel section.
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
const TemperatureIndex * model
VariableMetadata m_variable
void set_units(const std::string &units, const std::string &glaciological_units)
std::shared_ptr< TSDiagnostic > Ptr
Scalar diagnostic reporting a "flux".
void update_temp_mj(const IceModelVec2S &surfelev, const IceModelVec2S &lat, const IceModelVec2S &lon)
Updates mean July near-surface air temperature.
LocalMassBalance::DegreeDayFactors degree_day_factors(int i, int j, double latitude)
A PDD implementation which computes the local mass balance based on an expectation integral.
An alternative PDD implementation which simulates a random process to get the number of PDDs.
std::shared_ptr< atmosphere::AtmosphereModel > m_atmosphere
virtual DiagnosticList diagnostics_impl() const
virtual void write_model_state_impl(const File &output) const
The default (empty implementation).
virtual void define_model_state_impl(const File &output) const
The default (empty implementation).
static IceModelVec2S::Ptr allocate_melt(IceGrid::ConstPtr grid)
static IceModelVec2S::Ptr allocate_temperature(IceGrid::ConstPtr grid)
static IceModelVec2S::Ptr allocate_accumulation(IceGrid::ConstPtr grid)
virtual TSDiagnosticList ts_diagnostics_impl() const
static IceModelVec2S::Ptr allocate_runoff(IceGrid::ConstPtr grid)
const IceModelVec2S & accumulation() const
Returns accumulation.
virtual void init_impl(const Geometry &geometry)
The interface of PISM's surface models.
virtual const IceModelVec2S & runoff_impl() const
virtual const IceModelVec2S & accumulation_impl() const
IceModelVec2S m_firn_depth
firn depth
virtual TSDiagnosticList ts_diagnostics_impl() const
virtual void init_impl(const Geometry &geometry)
IceModelVec2S m_mass_flux
cached surface mass balance rate
virtual void update_impl(const Geometry &geometry, double t, double dt)
virtual const IceModelVec2S & mass_flux_impl() const
virtual const IceModelVec2S & melt_impl() const
TemperatureIndex(IceGrid::ConstPtr g, std::shared_ptr< atmosphere::AtmosphereModel > input)
std::shared_ptr< IceModelVec2T > m_air_temp_sd
standard deviation of the daily variability of the air temperature
std::unique_ptr< FaustoGrevePDDObject > m_faustogreve
if not NULL then user wanted fausto PDD stuff
double m_next_balance_year_start
virtual const IceModelVec2S & temperature_impl() const
double compute_next_balance_year_start(double time)
virtual MaxTimestep max_timestep_impl(double t) const
IceModelVec2S m_snow_depth
snow depth (reset once a year)
double m_base_pddStdDev
K; daily amount of randomness.
IceModelVec2S::Ptr m_temperature
const IceModelVec2S & firn_depth() const
virtual void define_model_state_impl(const File &output) const
The default (empty implementation).
IceModelVec2S::Ptr m_runoff
total runoff during the last time step
const IceModelVec2S & air_temp_sd() const
const IceModelVec2S & snow_depth() const
IceModelVec2S::Ptr m_accumulation
total accumulation during the last time step
virtual void write_model_state_impl(const File &output) const
The default (empty implementation).
virtual DiagnosticList diagnostics_impl() const
std::unique_ptr< LocalMassBalance > m_mbscheme
mass balance scheme to use
IceModelVec2S::Ptr m_melt
total melt during the last time step
LocalMassBalance::DegreeDayFactors m_base_ddf
holds degree-day factors in location-independent case
A class implementing a temperature-index (positive degree-day) scheme to compute melt and runoff,...
IceModelVec2S m_accumulation_mass
const IceModelVec2S & model_input()
Accumulation(const TemperatureIndex *m, AmountKind kind)
Report accumulation (precipitation minus rain), averaged over the reporting interval.
IceModelVec2S m_melt_mass
const IceModelVec2S & model_input()
SurfaceMelt(const TemperatureIndex *m, AmountKind kind)
Report surface melt, averaged over the reporting interval.
SurfaceRunoff(const TemperatureIndex *m, AmountKind kind)
IceModelVec2S m_runoff_mass
const IceModelVec2S & model_input()
Report surface runoff, averaged over the reporting interval.
TotalSurfaceAccumulation(const TemperatureIndex *m)
Reports the total accumulation rate.
TotalSurfaceMelt(const TemperatureIndex *m)
Reports the total melt rate.
TotalSurfaceRunoff(const TemperatureIndex *m)
Reports the total top surface ice flux.
#define PISM_ERROR_LOCATION
static double integrate(const IceModelVec2S &input)
double convert(System::Ptr system, double input, const std::string &spec1, const std::string &spec2)
Convert a quantity from unit1 to unit2.
InputOptions process_input_options(MPI_Comm com, Config::ConstPtr config)
std::map< std::string, TSDiagnostic::Ptr > TSDiagnosticList
std::map< std::string, Diagnostic::Ptr > DiagnosticList
@ 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)
double ice
m day^-1 K^-1; ice-equivalent amount of ice melted, per PDD
double refreeze_fraction
fraction of melted snow which refreezes as ice
double snow
m day^-1 K^-1; ice-equivalent amount of snow melted, per PDD
A struct which holds degree day factors.
static double S(unsigned n)