IceModel Class Reference

The base class for PISM. Contains all essential variables, parameters, and flags for modelling an ice sheet. More...

#include <iceModel.hh>

Inheritance diagram for IceModel:

Public Member Functions
	IceModel (IceGrid &g, NCConfigVariable &config, NCConfigVariable &overrides)
virtual	~IceModel ()
virtual PetscErrorCode	grid_setup ()
	Sets up the computational grid.
virtual PetscErrorCode	init_physics ()
	Initialize some physical parameters.
virtual PetscErrorCode	init_couplers ()
	Initializes atmosphere and ocean couplers.
virtual PetscErrorCode	set_grid_from_options ()
	Initalizes the grid from options.
virtual PetscErrorCode	set_grid_defaults ()
	Set default values of grid parameters.
virtual PetscErrorCode	model_state_setup ()
	Sets the starting values of model state variables.
virtual PetscErrorCode	set_vars_from_options ()
	Sets starting values of model state variables using command-line options.
virtual PetscErrorCode	report_grid_parameters ()
virtual PetscErrorCode	allocate_internal_objects ()
	Allocates work vectors (and calls more).
virtual PetscErrorCode	misc_setup ()
	Miscellaneous initialization tasks plus tasks that need the fields that can come from regridding.
PetscErrorCode	init ()
	Manage the initialization of the IceModel object.
virtual PetscErrorCode	run ()
	Do the time-stepping for an evolution run.
virtual PetscErrorCode	step (bool do_mass_continuity, bool do_energy, bool do_age, bool do_skip, bool use_ssa_when_grounded)
	The contents of the time-step.
virtual PetscErrorCode	setExecName (const char *my_executable_short_name)
virtual IceFlowLawFactory &	getIceFlowLawFactory ()
virtual IceFlowLaw *	getIceFlowLaw ()
virtual PetscErrorCode	bootstrapFromFile (const char *fname)
	Read file and use heuristics to initialize PISM from typical 2d data available through remote sensing.
virtual PetscErrorCode	readShelfStreamBCFromFile (const char *fname)
	Read certain boundary conditions from a NetCDF file, for diagnostic SSA calculations.
virtual PetscErrorCode	setFromOptions ()
	Read runtime (command line) options and alter the corresponding parameters or flags as appropriate.
virtual PetscErrorCode	set_output_size (string option, string description, string default_value, set< string > &result)
	Assembles a list of variables corresponding to an output file size.
virtual void	attach_surface_model (PISMSurfaceModel *surf)
virtual void	attach_ocean_model (PISMOceanModel *ocean)
virtual PetscErrorCode	additionalAtStartTimestep ()
	Virtual. Does nothing in `IceModel`. Derived classes can do more computation in each time step.
virtual PetscErrorCode	additionalAtEndTimestep ()
	Virtual. Does nothing in `IceModel`. Derived classes can do more computation in each time step.
virtual PetscErrorCode	correct_cell_areas ()
virtual PetscErrorCode	initFromFile (const char *)
	Read a saved PISM model state in NetCDF format, for complete initialization of an evolution or diagnostic run.
virtual PetscErrorCode	writeFiles (const char *default_filename)
	Save model state in NetCDF format.
virtual PetscErrorCode	write_model_state (const char *filename)
virtual PetscErrorCode	write_variables (const char *filename, set< string > vars, nc_type nctype)
	Writes variables listed in vars to filename, using nctype to write fields stored in dedicated IceModelVecs.
virtual PetscErrorCode	write_extra_fields (const char *filename)
	Writes extra fields to the output file `filename`. Does nothing in the base class.
Protected Member Functions
virtual PetscErrorCode	createVecs ()
	Allocate all IceModelVecs defined in IceModel.
virtual PetscErrorCode	deallocate_internal_objects ()
	De-allocate internal objects.
virtual void	setConstantNuHForSSA (PetscScalar)
virtual PetscErrorCode	computeMaxDiffusivity (bool update_diffusivity_viewer)
	Compute the maximum diffusivity associated to the SIA deformational velocity.
virtual PetscErrorCode	computeMax3DVelocities ()
	Compute the maximum velocities for time-stepping and reporting to user.
virtual PetscErrorCode	computeMax2DSlidingSpeed ()
	Compute the CFL constant associated to first-order upwinding for the sliding contribution to mass continuity.
virtual PetscErrorCode	adaptTimeStepDiffusivity ()
	Compute the maximum time step allowed by the diffusive SIA.
virtual PetscErrorCode	determineTimeStep (const bool doTemperatureCFL)
	Use various stability criteria to determine the time step for an evolution run.
virtual PetscErrorCode	countCFLViolations (PetscScalar *CFLviol)
	Because of the -skip mechanism it is still possible that we can have CFL violations: count them.
virtual PetscErrorCode	initBasalTillModel ()
	Initialize the pseudo-plastic till mechanical model.
virtual PetscErrorCode	computePhiFromBedElevation ()
	Computes the till friction angle phi as a piecewise linear function of bed elevation, according to user options.
virtual PetscScalar	getBasalWaterPressure (PetscScalar thk, PetscScalar bwat, PetscScalar bmr, PetscScalar frac, PetscScalar hmelt_max) const
	Compute modeled pressure in subglacial liquid water using thickness of water layer, and possibly the melt rate, at the base.
virtual PetscErrorCode	updateYieldStressUsingBasalWater ()
	Update the till yield stress for the pseudo-plastic till SSA model.
virtual PetscScalar	basalDragx (PetscScalar tauc, PISMVector2 uv, PetscInt i, PetscInt j) const
virtual PetscScalar	basalDragy (PetscScalar tauc, PISMVector2 uv, PetscInt i, PetscInt j) const
virtual PetscErrorCode	diffuseHmelt ()
	Apply explicit time step for pure diffusion to basal layer of melt water.
virtual PetscErrorCode	bed_def_setup ()
virtual PetscErrorCode	bed_def_step ()
virtual PetscErrorCode	putTempAtDepth ()
	Create a temperature field within ice and bedrock from given surface temperature and geothermal flux maps.
virtual PetscErrorCode	bootstrapSetBedrockColumnTemp (PetscInt i, PetscInt j, PetscScalar Ttopbedrock, PetscScalar geothermflux, PetscScalar bed_thermal_k)
	Set the temperatures in a column of bedrock based on a temperature at the top and a geothermal flux.
virtual PetscErrorCode	setMaskSurfaceElevation_bootstrap ()
	Determine surface and mask according to information in bootstrap file and options.
PetscErrorCode	setDefaults ()
	Assigns default values to the many parameters and flags in IceModel.
virtual PetscErrorCode	setEnth3FromT3_ColdIce ()
	Compute Enth3 from temperature T3 by assuming the ice has zero liquid fraction.
virtual PetscErrorCode	setEnth3FromT3AndLiqfrac3 (IceModelVec3 &Liqfrac3)
	Compute Enth3 from temperature T3 and liquid fraction.
virtual PetscErrorCode	setLiquidFracFromEnthalpy (IceModelVec3 &useForLiquidFrac)
	Compute the liquid fraction corresponding to Enth3, and put in a global IceModelVec3 provided by user.
virtual PetscErrorCode	setCTSFromEnthalpy (IceModelVec3 &useForCTS)
	Compute the CTS field, CTS = E/E_s(p), from Enth3, and put in a global IceModelVec3 provided by user.
virtual PetscErrorCode	getEnthalpyCTSColumn (PetscScalar p_air, PetscScalar thk, PetscInt ks, const PetscScalar Enth, const PetscScalar w, PetscScalar lambda, PetscScalar *Enth_s)
	Compute the CTS value of enthalpy in an ice column, and the lambda for BOMBPROOF.
virtual PetscErrorCode	enthalpyAndDrainageStep (PetscScalar vertSacrCount, PetscScalar liquifiedVol)
	Update enthalpy field based on conservation of energy in ice and bedrock.
virtual PetscErrorCode	computeDrivingStress (IceModelVec2S &vtaudx, IceModelVec2S &vtaudy)
	Compute vector driving stress at base of ice on the regular grid.
virtual PetscErrorCode	updateSurfaceElevationAndMask ()
	Update the surface elevation and the flow-type mask when the geometry has changed.
virtual PetscErrorCode	update_mask ()
virtual PetscErrorCode	update_surface_elevation ()
virtual PetscErrorCode	massContExplicitStep ()
	Update the thickness from the horizontal velocity and the surface and basal mass balance.
virtual PetscScalar	grainSizeVostok (PetscScalar age) const
	Use the Vostok core as a source of a relationship between the age of the ice and the grain size.
virtual PetscErrorCode	set_time_from_options ()
	Determine the run length, starting and ending years using command-line options.
virtual PetscErrorCode	dumpToFile (const char *filename)
virtual PetscErrorCode	regrid ()
	Manage regridding based on user options. Call IceModelVec.regrid() to do each selected variable.
virtual PetscErrorCode	writeSSAsystemMatlab (IceModelVec2S vNuH[2])
	Write out the linear system of equations for the SSA linear solve at the last nonlinear iteration in a time step.
virtual PetscErrorCode	computeFlowUbarStats (PetscScalar gUbarmax, PetscScalar gUbarSIAav, PetscScalar gUbarstreamav, PetscScalar gUbarshelfav, PetscScalar gicegridfrac, PetscScalar gSIAgridfrac, PetscScalar gstreamgridfrac, PetscScalar gshelfgridfrac)
virtual PetscErrorCode	volumeArea (PetscScalar &gvolume, PetscScalar &garea, PetscScalar &gvolSIA, PetscScalar &gvolstream, PetscScalar &gvolshelf)
	Computes volume and area of ice sheet, for reporting purposes.
virtual PetscErrorCode	energyStats (PetscScalar iarea, bool useHomoTemp, PetscScalar &gmeltfrac, PetscScalar &gtemp0)
virtual PetscErrorCode	ageStats (PetscScalar ivol, PetscScalar &gorigfrac)
virtual PetscErrorCode	summary (bool tempAndAge, bool useHomoTemp)
virtual PetscErrorCode	summaryPrintLine (PetscTruth printPrototype, bool tempAndAge, PetscScalar year, PetscScalar delta_t, PetscScalar volume, PetscScalar area, PetscScalar meltfrac, PetscScalar H0, PetscScalar T0)
	Print a line to stdout which summarizes the state of the modeled ice sheet at the end of the time step.
virtual PetscErrorCode	compute_by_name (string name, IceModelVec *&result)
	Computes a diagnostic quantity given by `name` and returns a pointer to a pre-allocated work vector containing it.
virtual PetscErrorCode	compute_bwp (IceModelVec2S &result)
	Computes the subglacial (basal) water pressure.
virtual PetscErrorCode	compute_cbar (IceModelVec2S &result)
	Computes cbar, the magnitude of vertically-integrated horizontal velocity of ice and masks out ice-free areas.
virtual PetscErrorCode	compute_cbase (IceModelVec2S &result, IceModelVec2S &tmp)
	Computes cbase, the magnitude of horizontal velocity of ice at base of ice and masks out ice-free areas. Uses `tmp` as a preallocated temporary storage.
virtual PetscErrorCode	compute_cflx (IceModelVec2S &result, IceModelVec2S &cbar)
	Computes cflx, the magnitude of vertically-integrated horizontal flux of ice.
virtual PetscErrorCode	compute_csurf (IceModelVec2S &result, IceModelVec2S &tmp)
	Computes csurf, the magnitude of horizontal velocity of ice at ice surface and masks out ice-free areas. Uses `tmp` as a preallocated temporary storage.
virtual PetscErrorCode	compute_dhdt (IceModelVec2S &result)
	Computes the rate of change of ice surface elevation as a sum of the bedrock uplift rate and the thickness rate of change.
virtual PetscErrorCode	compute_enthalpybase (IceModelVec2S &result)
	Computes ice enthalpy at the base of ice.
virtual PetscErrorCode	compute_enthalpysurf (IceModelVec2S &result)
	Computes ice enthalpy at the 1 m below the surface.
virtual PetscErrorCode	compute_hardav (IceModelVec2S &result)
	Computes the vertically-averaged ice hardness.
virtual PetscErrorCode	compute_taud (IceModelVec2S &result, IceModelVec2S &tmp)
	Computes taud, the magnitude of driving shear stress at base of ice. Uses tmp as a preallocated temporary storage.
virtual PetscErrorCode	compute_cts (IceModelVec3 &useForCTS)
	Compute the CTS field, CTS = E/E_s(p) and put in a global IceModelVec3 provided by user.
virtual PetscErrorCode	compute_liqfrac (IceModelVec3 &useForLiqfrac)
	Compute the liquid fraction, and put in a global IceModelVec3 provided by user.
virtual PetscErrorCode	compute_temp (IceModelVec3 &result)
	Compute the ice temperature corresponding to Enth3, and put in a given IceModelVec3.
virtual PetscErrorCode	compute_temp_pa (IceModelVec3 &result)
	Compute the pressure-adjusted temperature in degrees C corresponding to ice temperature, and put in a global IceModelVec3 provided by user.
virtual PetscErrorCode	compute_tempbase (IceModelVec2S &result)
	Computes ice temperature at the base of ice.
virtual PetscErrorCode	compute_temppabase (IceModelVec3 &hasPATemp, IceModelVec2S &result)
	Computes pressure-adjusted ice temperature at the base of ice.
virtual PetscErrorCode	compute_tempsurf (IceModelVec2S &result)
	Computes ice temperature at the 1 m below the surface.
virtual PetscErrorCode	compute_uvelbase (IceModelVec2S &result)
virtual PetscErrorCode	compute_uvelsurf (IceModelVec2S &result)
	Computes uvelsurf, the x component of velocity of ice at ice surface.
virtual PetscErrorCode	compute_vvelbase (IceModelVec2S &result)
virtual PetscErrorCode	compute_vvelsurf (IceModelVec2S &result)
	Computes vvelsurf, the y component of velocity of ice at ice surface.
virtual PetscErrorCode	compute_wvelbase (IceModelVec2S &result)
virtual PetscErrorCode	compute_wvelsurf (IceModelVec2S &result)
	Computes wvelsurf, the vertical velocity of ice at ice surface.
virtual PetscErrorCode	compute_rank (IceModelVec2S &result)
	Sets entrues of result to corresponding processor ranks.
virtual PetscErrorCode	compute_proc_ice_area (IceModelVec2S &result)
	Sets entrues of result to corresponding processor ranks.
virtual PetscErrorCode	ice_mass_bookkeeping ()
virtual PetscErrorCode	compute_ice_volume (PetscScalar &result)
	Computes the ice volume, in m^3.
virtual PetscErrorCode	compute_ice_area (PetscScalar &result)
	Computes ice area, in m^2.
virtual PetscErrorCode	compute_ice_area_grounded (PetscScalar &result)
	Computes grounded ice area, in m^2.
virtual PetscErrorCode	compute_ice_area_floating (PetscScalar &result)
	Computes floating ice area, in m^2.
virtual PetscErrorCode	compute_ice_enthalpy (PetscScalar &result)
	Computes the total ice enthalpy in J.
virtual PetscErrorCode	compute_by_name (string name, PetscScalar &result)
	Compute a scalar diagnostic quantity by name.
virtual PetscErrorCode	surfaceGradientSIA ()
	Compute the surface gradient in advance of the SIA velocity computation.
virtual PetscErrorCode	velocitySIAStaggered ()
	Compute the vertically-averaged horizontal velocity according to the non-sliding SIA.
virtual PetscErrorCode	basalSlidingHeatingSIA ()
	Compute the basal sliding and frictional heating if (where) SIA sliding rule is used.
virtual PetscErrorCode	velocities2DSIAToRegular ()
	Average staggered-grid vertically-averaged horizontal velocity onto regular grid.
virtual PetscErrorCode	SigmaSIAToRegular ()
	Put the volume strain heating (dissipation heating) onto the regular grid.
virtual PetscErrorCode	horizontalVelocitySIARegular ()
	Update regular grid horizontal velocities u,v at depth for SIA regions.
virtual PetscScalar	basalVelocitySIA (PetscScalar x, PetscScalar y, PetscScalar H, PetscScalar T, PetscScalar alpha, PetscScalar mu, PetscScalar min_T) const
	Compute the coefficient of surface gradient, for basal sliding velocity as a function of driving stress in SIA regions.
virtual PetscErrorCode	allocateSSAobjects ()
	Allocates SSA tools.
virtual PetscErrorCode	destroySSAobjects ()
	Deallocate SSA tools.
virtual PetscErrorCode	initSSA ()
	Each step of SSA uses previously saved values to start iteration; zero them here to start.
virtual PetscErrorCode	velocitySSA (PetscInt *numiter)
	Compute the vertically-averaged horizontal velocity from the shallow shelf approximation (SSA).
virtual PetscErrorCode	velocitySSA (IceModelVec2S vNuH[2], PetscInt *numiter)
	Call this one directly if control over allocation of vNuH[2] is needed (e.g. test J).
virtual PetscErrorCode	computeEffectiveViscosity (IceModelVec2S vNuH[2], PetscReal epsilon)
	Compute the product of the effective viscosity $\nu$ and ice thickness for the SSA model.
virtual PetscErrorCode	testConvergenceOfNu (IceModelVec2S vNuH[2], IceModelVec2S vNuHOld[2], PetscReal norm, PetscReal normChange)
virtual PetscErrorCode	assembleSSAMatrix (bool includeBasalShear, IceModelVec2S vNuH[2], Mat A)
	Assemble the left-hand side matrix for the numerical approximation of the SSA equations.
virtual PetscErrorCode	assembleSSARhs (Vec rhs)
	Computes the right-hand side ("rhs") of the linear problem for the SSA equations.
virtual PetscErrorCode	trivialMoveSSAXtoIMV2V ()
virtual PetscErrorCode	broadcastSSAVelocity (bool updateVelocityAtDepth)
	At all SSA points, update the velocity field.
virtual PetscErrorCode	correctSigma ()
	At SSA points, correct the previously-computed volume strain heating (dissipation heating).
virtual PetscErrorCode	correctBasalFrictionalHeating ()
	At SSA points, correct the previously-computed basal frictional heating.
virtual PetscErrorCode	energyStep ()
	Manage the solution of the energy equation, and related parallel communication.
virtual PetscErrorCode	temperatureStep (PetscScalar vertSacrCount, PetscScalar bulgeCount)
	Takes a semi-implicit time-step for the temperature equation.
virtual PetscErrorCode	ageStep ()
	Take a semi-implicit time-step for the age equation.
virtual bool	checkThinNeigh (PetscScalar E, PetscScalar NE, PetscScalar N, PetscScalar NW, PetscScalar W, PetscScalar SW, PetscScalar S, PetscScalar SE)
virtual PetscErrorCode	excessToFromBasalMeltLayer (const PetscScalar rho_c, const PetscScalar z, const PetscScalar dz, PetscScalar Texcess, PetscScalar Hmelt)
	Compute the melt water which should go to the base if is above pressure-melting.
virtual int	endOfTimeStepHook ()
	Catch signals -USR1 and -TERM; in the former case save and continue; in the latter, save and stop.
virtual PetscErrorCode	stampHistoryCommand ()
	Build a history string from the command which invoked PISM.
virtual PetscErrorCode	stampHistoryEnd ()
	Build the particular history string associated to the end of a PISM run.
virtual PetscErrorCode	stampHistory (string)
	Get time and user/host name and add it to the given string.
virtual PetscErrorCode	check_maximum_thickness ()
	Check if the thickness of the ice is too large and extend the grid if necessary.
virtual PetscErrorCode	check_maximum_thickness_hook (const int old_Mz)
	Allows derived classes to extend their own IceModelVec3's in vertical.
virtual bool	issounding (const PetscInt i, const PetscInt j)
virtual PetscErrorCode	velocity (bool updateSIAVelocityAtDepth)
	Manage the computation of velocity and do the necessary parallel communication.
virtual PetscErrorCode	vertVelocityFromIncompressibility ()
	Compute vertical velocity using incompressibility of the ice.
PetscErrorCode	init_snapshots ()
	Initializes the snapshot-saving mechanism.
PetscErrorCode	write_snapshot ()
	Writes a snapshot of the model state (if necessary).
PetscErrorCode	init_timeseries ()
	Initializes the code writing scalar time-series.
PetscErrorCode	create_timeseries ()
	Creates DiagnosticTimeseries objects used to store and report scalar diagnostic quantities.
PetscErrorCode	write_timeseries ()
	Write time-series.
PetscErrorCode	ts_max_timestep (double t_years, double &dt_years)
	Computes the maximum time-step we can take and still hit all the requested years.
PetscErrorCode	init_extras ()
	Initialize the code saving spatially-variable diagnostic quantities.
PetscErrorCode	write_extras ()
	Write spatially-variable diagnostic quantities.
PetscErrorCode	extras_max_timestep (double t_years, double &dt_years)
	Computes the maximum time-step we can take and still hit all the requested years.
virtual PetscErrorCode	init_viewers ()
	Initialize run-time diagnostic viewers.
virtual PetscErrorCode	update_viewers ()
	Update the runtime graphical viewers.
virtual PetscErrorCode	update_nu_viewers (IceModelVec2S vNu[2], IceModelVec2S[2], bool)
	Update nuH viewers.
Protected Attributes
IceGrid &	grid
NCConfigVariable	mapping
	grid projection (mapping) parameters
NCConfigVariable &	config
	configuration flags and parameters
NCConfigVariable &	overrides
	flags and parameters overriding config, see -config_override
NCGlobalAttributes	global_attributes
IceFlowLawFactory	iceFactory
IceFlowLaw *	ice
IceBasalResistancePlasticLaw *	basal
SSAStrengthExtension	ssaStrengthExtend
EnthalpyConverter *	EC
PISMSurfaceModel *	surface
PISMOceanModel *	ocean
PISMBedDef *	beddef
PISMVars	variables
	A dictionary with pointers to IceModelVecs below, for passing them from the IceModel core to other components (such as surface and ocean models).
IceModelVec2S	vh
	ice surface elevation; ghosted
IceModelVec2S	vH
	ice thickness; ghosted
IceModelVec2S	vdHdt
	$\frac{dH}{dt}$ ; ghosted to simplify the code computing it
IceModelVec2S	vtauc
	yield stress for basal till (plastic or pseudo-plastic model); no ghosts
IceModelVec2S	vHmelt
	thickness of the basal meltwater; ghosted (because of the diffusion)
IceModelVec2S	vbmr
	rate of production of basal meltwater (ice-equivalent); no ghosts
IceModelVec2S	vLongitude
	Longitude; ghosted to compute cell areas.
IceModelVec2S	vLatitude
	Latitude; ghosted to compute cell areas.
IceModelVec2S	vbed
	bed topography; ghosted
IceModelVec2S	vuplift
	bed uplift rate; ghosted to simplify the code computing it
IceModelVec2S	vGhf
	geothermal flux; no ghosts
IceModelVec2S	vRb
	basal frictional heating on regular grid; no ghosts
IceModelVec2S	vtillphi
	friction angle for till under grounded ice sheet; no ghosts
IceModelVec2S	acab
	accumulation/ablation rate; no ghosts
IceModelVec2S	artm
	ice temperature at the ice surface but below firn; no ghosts
IceModelVec2S	shelfbtemp
	ice temperature at the shelf base; no ghosts
IceModelVec2S	shelfbmassflux
	ice mass flux into the ocean at the shelf base; no ghosts
IceModelVec2S	cell_area
	cell areas (computed using the WGS84 datum)
IceModelVec2Stag	uvbar
	ubar and vbar on staggered grid; ubar at i+1/2, vbar at j+1/2
IceModelVec2V	vel_basal
	basal velocities on standard grid; ghosted
IceModelVec2V	vel_bar
	vertically-averaged horizontal velocity on standard grid; ghosted
IceModelVec2Mask	vMask
	mask for flow type with values SHEET, DRAGGING, FLOATING
IceModelVec3	u3
IceModelVec3	v3
IceModelVec3	w3
	velocity of ice; m s-1 (ghosted)
IceModelVec3	Sigma3
	strain-heating term in conservation of energy model; J s-1 m-3 (no ghosts)
IceModelVec3	T3
	absolute temperature of ice; K (ghosted)
IceModelVec3	Enth3
	enthalpy; J / kg (ghosted)
IceModelVec3	tau3
	age of ice; s (ghosted because it is evaluated on the staggered-grid)
IceModelVec3Bedrock	Tb3
	temperature of lithosphere (bedrock) under ice or ocean; K (no ghosts)
PetscReal	dt
	mass continuity time step, s
PetscReal	dtTempAge
PetscReal	maxdt_temporary
PetscReal	dt_force
PetscReal	CFLviolcount
	really is just a count, but PetscGlobalSum requires this type
PetscReal	dt_from_diffus
PetscReal	dt_from_cfl
PetscReal	CFLmaxdt
PetscReal	CFLmaxdt2D
PetscReal	gDmax
PetscReal	gmaxu
PetscReal	gmaxv
PetscReal	gmaxw
PetscReal	gdHdtav
	average value in map-plane (2D) of dH/dt, where there is ice; m s-1
PetscReal	total_sub_shelf_ice_flux
PetscReal	total_basal_ice_flux
PetscReal	total_surface_ice_flux
PetscReal	dvoldt
	d(total ice volume)/dt; m3 s-1
PetscInt	skipCountDown
PetscScalar	standard_gravity
bool	leaveNuHAloneSSA
PetscTruth	updateHmelt
PetscTruth	holdTillYieldStress
PetscTruth	useConstantTillPhi
PetscTruth	shelvesDragToo
PetscTruth	doAdaptTimeStep
PetscTruth	realAgeForGrainSize
PetscTruth	ssaSystemToASCIIMatlab
PetscTruth	reportPATemps
PetscTruth	allowAboveMelting
PetscTruth	computeSIAVelocities
char	adaptReasonFlag
char	ssaMatlabFilePrefix [PETSC_MAX_PATH_LEN]
string	stdout_flags
string	stdout_ssa
string	executable_short_name
PetscScalar	last_bed_def_update
IceModelVec2S	vWork2d [nWork2d]
IceModelVec3	Tnew3
IceModelVec3	Enthnew3
IceModelVec3	taunew3
IceModelVec3	Sigmastag3 [2]
IceModelVec3	Istag3 [2]
int	have_ssa_velocities
IceModelVec2V	vel_ssa
IceModelVec2V	vel_ssa_old
KSP	SSAKSP
Mat	SSAStiffnessMatrix
Vec	SSAX
Vec	SSARHS
DA	SSADA
set< string >	output_vars
string	snapshots_filename
bool	save_snapshots
bool	snapshots_file_is_ready
bool	split_snapshots
vector< double >	snapshot_times
set< string >	snapshot_vars
unsigned int	current_snapshot
bool	save_ts
string	ts_filename
	true if the user requested time-series output
vector< double >	ts_times
	file to write time-series to
unsigned int	current_ts
	times requested
set< string >	ts_vars
	index of the current time
vector< DiagnosticTimeseries * >	timeseries
	variables requested
bool	save_extra
bool	extra_file_is_ready
bool	split_extra
string	extra_filename
vector< double >	extra_times
unsigned int	current_extra
set< string >	extra_vars
set< string >	map_viewers
set< string >	slice_viewers
set< string >	sounding_viewers
PetscInt	id
PetscInt	jd
bool	view_diffusivity
bool	view_log_nuH
bool	view_nuH
Static Protected Attributes
static const PetscInt	nWork2d = 5
Private Attributes
int	siaEVENT
int	ssaEVENT
int	velmiscEVENT
int	beddefEVENT
int	massbalEVENT
int	tempEVENT

Detailed Description

The base class for PISM. Contains all essential variables, parameters, and flags for modelling an ice sheet.

Definition at line 55 of file iceModel.hh.

Constructor & Destructor Documentation

IceModel	(	IceGrid &	g,
		NCConfigVariable &	config,
		NCConfigVariable &	overrides
	)

Definition at line 26 of file iceModel.cc.

References allowAboveMelting, basal, beddef, CFLviolcount, IceGrid.com, doAdaptTimeStep, dvoldt, EC, gdHdtav, global_attributes, grid, have_ssa_velocities, ICE_PB, iceFactory, NCVariable.init(), mapping, ocean, pism_signal, pism_signal_handler(), IceGrid.rank, save_extra, save_snapshots, save_ts, setDefaults(), IceFlowLawFactory.setType(), surface, total_basal_ice_flux, total_sub_shelf_ice_flux, total_surface_ice_flux, verbPrintf(), view_diffusivity, view_log_nuH, and view_nuH.

~IceModel ( ) [virtual]

Definition at line 86 of file iceModel.cc.

References basal, beddef, deallocate_internal_objects(), EC, ice, ocean, surface, and timeseries.

Member Function Documentation

PetscErrorCode adaptTimeStepDiffusivity ( ) [protected, virtual]

Compute the maximum time step allowed by the diffusive SIA.

Note computeMaxDiffusivity() must be called before this to set gDmax. Note adapt_ratio * 2 is multiplied by dx^2/(2*maxD) so dt <= adapt_ratio * dx^2/maxD (if dx=dy)

Reference: MortonMayers pp 62--63.

Definition at line 240 of file iMadaptive.cc.

References adaptReasonFlag, config, dt, dt_from_diffus, IceGrid.dx, IceGrid.dy, gDmax, NCConfigVariable.get(), NCConfigVariable.get_flag(), grid, and skipCountDown.

Referenced by determineTimeStep().

PetscErrorCode additionalAtEndTimestep ( ) [virtual]

Virtual. Does nothing in IceModel. Derived classes can do more computation in each time step.

Reimplemented in IcePSTexModel, IceMISMIPModel, and IceCompModel.

Definition at line 32 of file iMutil.cc.

Referenced by step().

PetscErrorCode additionalAtStartTimestep ( ) [virtual]

Virtual. Does nothing in IceModel. Derived classes can do more computation in each time step.

Reimplemented in IceMISMIPModel, and IceCompModel.

Definition at line 26 of file iMutil.cc.

Referenced by step().

PetscErrorCode ageStats	(	PetscScalar	ivol,
		PetscScalar &	gorigfrac
	)			`[protected, virtual]`

Computes fraction of the ice which is as old as the start of the run (original). Communication occurs here.

Definition at line 203 of file iMreport.cc.

References IceModelVec.begin_access(), IceGrid.com, config, IceGrid.dx, IceGrid.dy, e, IceModelVec.end_access(), IceModelVec2S.get_array(), NCConfigVariable.get_flag(), IceModelVec3.getInternalColumn(), grid, IceGrid.kBelowHeight(), secpera, tau3, vH, IceGrid.xm, IceGrid.xs, IceGrid.year, IceGrid.ym, IceGrid.ys, and IceGrid.zlevels.

Referenced by summary().

PetscErrorCode ageStep ( ) [protected, virtual]

Take a semi-implicit time-step for the age equation.

The age equation is $d\tau/dt = 1$ , that is,

$\frac{\partial \tau}{\partial t} + u \frac{\partial \tau}{\partial x} + v \frac{\partial \tau}{\partial y} + w \frac{\partial \tau}{\partial z} = 1$

where $\tau(t,x,y,z)$ is the age of the ice and $(u,v,w)$ is the three dimensional velocity field. This equation is purely advective. And it is hyperbolic.

The boundary condition is that when the ice falls as snow it has age zero. That is, $\tau(t,x,y,h(t,x,y)) = 0$ in accumulation areas, while there is no boundary condition elsewhere, as the characteristics go outward in the ablation zone. (Some more numerical-analytic attention to this is worthwhile.)

If the velocity in the bottom cell of ice is upward (

 (w[i][j][0] > 0

) then we also apply an age = 0 boundary condition. This is the case where ice freezes on at the base, either grounded basal ice freezing on stored water in till, or marine basal ice.

The numerical method is first-order upwind but the vertical advection term is computed implicitly. (Thus there is no CFL-type stability condition for that part.)

We use a finely-spaced, equally-spaced vertical grid in the calculation. Note that the IceModelVec3 methods getValColumn...() and setValColumn..() interpolate back and forth between the grid on which calculation is done and the storage grid. Thus the storage grid can be either equally spaced or not.

Definition at line 565 of file iMtemp.cc.

References IceModelVec.begin_access(), IceModelVec3.beginGhostCommTransfer(), IceGrid.com, columnSystemCtx.ddratio(), ageSystemCtx.dtAge, dtTempAge, IceGrid.dx, ageSystemCtx.dx, IceGrid.dy, ageSystemCtx.dy, IceGrid.dz_fine, ageSystemCtx.dzEQ, IceModelVec.end_access(), IceModelVec3.endGhostCommTransfer(), IceModelVec2S.get_array(), IceModelVec3.getValColumn(), grid, ageSystemCtx.initAllColumns(), IceGrid.Mz_fine, columnSystemCtx.norm1(), IceModelVec3.setColumn(), columnSystemCtx.setIndicesAndClearThisColumn(), IceModelVec3.setValColumnPL(), ageSystemCtx.solveThisColumn(), tau3, ageSystemCtx.tau3, taunew3, ageSystemCtx.u, u3, ageSystemCtx.v, v3, vH, ageSystemCtx.w, w3, x, IceGrid.xm, IceGrid.xs, IceGrid.ym, IceGrid.ys, and IceGrid.zlevels_fine.

Referenced by step().

PetscErrorCode allocate_internal_objects ( ) [virtual]

Allocates work vectors (and calls more).

Definition at line 484 of file iMinit.cc.

References allocateSSAobjects(), config, IceModelVec3.create(), IceModelVec2V.create(), IceModelVec2S.create(), Enthnew3, NCConfigVariable.get_flag(), grid, Istag3, MY_WIDE_STENCIL, nWork2d, IceModelVec.set_attrs(), Sigmastag3, taunew3, Tnew3, vel_ssa_old, and vWork2d.

Referenced by init().

PetscErrorCode allocateSSAobjects ( ) [protected, virtual]

Allocates SSA tools.

Definition at line 30 of file iMssa.cc.

References IceGrid.com, grid, IceGrid.Mx, IceGrid.My, IceGrid.Nx, IceGrid.Ny, SSADA, SSAKSP, SSARHS, SSAStiffnessMatrix, and SSAX.

Referenced by allocate_internal_objects().

PetscErrorCode assembleSSAMatrix	(	bool	includeBasalShear,
		IceModelVec2S	vNuH[2],
		Mat	A
	)			`[protected, virtual]`

Assemble the left-hand side matrix for the numerical approximation of the SSA equations.

The SSA equations are in their clearest divergence form

$- \frac{\partial T_{ij}}{\partial x_j} - \tau_{(b)i} = f_i$

where $i,j$ range over $x,y$ , $T_{ij}$ is a depth-integrated viscous stress tensor (i.e. equation (2.6) in [SchoofStream]; also [Morland]). These equations determine velocity in a more-or-less elliptic manner. Here $\tau_{(b)i}$ are the components of the basal shear stress applied to the base of the ice and $f_i$ is the driving shear stress,

$f_i = - \rho g H \frac{\partial h}{\partial x_i}.$

Here $H$ is the ice thickness and $h$ is the elevation of the surface of the ice.

More concretely, the SSA equations are

$\begin{align*} - 2 \left[\bar\nu H \left(2 u_x + v_y\right)\right]_x - \left[\bar\nu H \left(u_y + v_x\right)\right]_y - \tau_{(b)1} &= - \rho g H h_x, \\ - \left[\bar\nu H \left(u_y + v_x\right)\right]_x - 2 \left[\bar\nu H \left(u_x + 2 v_y\right)\right]_y - \tau_{(b)2} &= - \rho g H h_y, \end{align*}$

where $u$ is the $x$ -component of the velocity and $v$ is the $y$ -component of the velocity. Note $\bar\nu$ is the vertically-averaged effective viscosity of the ice.

For ice shelves $\tau_{(b)i} = 0$ [MacAyealetal]. For ice streams with a basal till modelled as a plastic material, $\tau_{(b)i} = - \tau_c u_i/|\mathbf{u}|$ where $\mathbf{u} = (u,v)$ , $|\mathbf{u}| = \left(u^2 + v^2\right)^{1/2}$ . Here $\tau_c(t,x,y)$ is the yield stress of the till [SchoofStream]. More generally, ice streams can be modeled with a pseudo-plastic basal till; see initBasalTillModel() and updateYieldStressUsingBasalWater() and [BKAJS].

The pseudo-plastic till model includes all power law sliding relations [BKAJS], and in particular it includes modeling the basal till as a linearly-viscous material, $\tau_{(b)i} = - \beta u_i$ where $\beta$ is the basal drag (friction) parameter [MacAyeal]. PISM assumes that the basal shear stress can be factored this way, even if the coefficient depends on the velocity, $\beta(u,v)$ . Such factoring is possible even in the case of (regularized) plastic till. This scalar coefficient $\beta$ is what is returned by IceBasalResistancePlasticLaw.drag().

Note that the basal shear stress appears on the left side of the linear system we actually solve. We believe this is crucial, because of its effect on the spectrum of the linear approximations of each stage. The effect on spectrum is clearest in the linearly-viscous till case but there seems to be an analogous effect in the plastic till case.

This method assembles the matrix for the left side of the above SSA equations. The numerical method is finite difference. Suppose we use difference notation $\delta_{+x}f^{i,j} = f^{i+1,j}-f^{i,j}$ , $\delta_{-x}f^{i,j} = f^{i,j}-f^{i-1,j}$ , and $\Delta_{x}f^{i,j} = f^{i+1,j}-f^{i-1,j}$ , and corresponding notation for $y$ differences, and that we write $N = \bar\nu$ then the first of the two "concrete" SSA equations above has this discretization:

$\begin{align*} - &2 \frac{N^{i+\frac{1}{2},j}}{\Delta x} \left[2\frac{\delta_{+x}u^{i,j}}{\Delta x} + \frac{\Delta_{y} v^{i+1,j} + \Delta_{y} v^{i,j}}{4 \Delta y}\right] + 2 \frac{N^{i-\frac{1}{2},j}}{\Delta x} \left[2\frac{\delta_{-x}u^{i,j}}{\Delta x} + \frac{\Delta_y v^{i,j} + \Delta_y v^{i-1,j}}{4 \Delta y}\right] \\ &\qquad- \frac{N^{i,j+\frac{1}{2}}}{\Delta y} \left[\frac{\delta_{+y} u^{i,j}}{\Delta y} + \frac{\Delta_x v^{i,j+1} + \Delta_x v^{i,j}}{4 \Delta x}\right] + \frac{N^{i,j-\frac{1}{2}}}{\Delta y} \left[\frac{\delta_{-y}u^{i,j}}{\Delta y} + \frac{\Delta_x v^{i,j} + \Delta_x v^{i,j-1}}{4 \Delta x}\right] - \tau_{(b)1}^{i,j} = - \rho g H^{i,j} \frac{\Delta_x h^{i,j}}{2\Delta x}. \end{align*}$

As a picture, see Figure ssastencil.

ssastencil: Stencil for our finite difference discretization of the first of the two scalar SSA equations. Triangles show staggered grid points where N = nu * H is evaluated. Circles and squares show where u and v are approximated, respectively.

It follows immediately that the matrix we assemble in the current method has 13 nonzeros entries per row because, for this first SSA equation, there are 5 grid values of $u$ and 8 grid values of $v$ used in this scheme. For the second equation we also have 13 nonzeros per row.

FIXME: address use of DAGetMatrix and MatStencil once those are used

Definition at line 385 of file iMssa.cc.

References basalDragx(), basalDragy(), IceModelVec.begin_access(), config, IceGrid.dx, IceGrid.dy, IceModelVec.end_access(), NCConfigVariable.get(), IceModelVec2V.get_array(), IceModelVec2S.get_array(), grid, MASK_DRAGGING_SHEET, MASK_FLOATING, MASK_SHEET, shelvesDragToo, IceModelVec2Mask.value(), vel_ssa, vMask, vtauc, IceGrid.xm, IceGrid.xs, IceGrid.ym, and IceGrid.ys.

Referenced by velocitySSA().

PetscErrorCode assembleSSARhs ( Vec rhs ) [protected, virtual]

Computes the right-hand side ("rhs") of the linear problem for the SSA equations.

The right side of the SSA equations is just the driving stress term

$- \rho g H \nabla h.$

The basal stress is put on the left side of the system. This method builds the discrete approximation of the right side. For more about the discretization of the SSA equations, see comments for assembleSSAMatrix().

The values of the driving stress on the i,j grid come from a call to computeDrivingStress().

Grid points with mask value MASK_SHEET correspond to the trivial equations

$\bar u_{ij} = \frac{uvbar(i-1,j,0) + uvbar(i,j,0)}{2},$

and similarly for $\bar v_{ij}$ . That is, the vertically-averaged horizontal velocity is already known for these points because it was either computed (on the staggered grid) using the SIA or was set by the -ssaBC mechanism.

Definition at line 561 of file iMssa.cc.

References IceModelVec2Stag.begin_access(), IceModelVec.begin_access(), computeDrivingStress(), IceModelVec2Stag.end_access(), IceModelVec.end_access(), IceModelVec2S.get_array(), grid, MASK_SHEET, SSADA, PISMVector2.u, uvbar, PISMVector2.v, IceModelVec2Mask.value(), vMask, vWork2d, IceGrid.xm, IceGrid.xs, IceGrid.ym, and IceGrid.ys.

Referenced by velocitySSA().

void attach_ocean_model ( PISMOceanModel * ocean ) [virtual]

Definition at line 365 of file iMutil.cc.

References ocean.

Referenced by main().

void attach_surface_model ( PISMSurfaceModel * surf ) [virtual]

Definition at line 361 of file iMutil.cc.

References surface.

Referenced by main().

PetscScalar basalDragx	(	PetscScalar **	tauc,
		PISMVector2 **	uv,
		PetscInt	i,
		PetscInt	j
	)			const `[protected, virtual]`

Reimplemented in IceMISMIPModel.

Definition at line 30 of file iMbasal.cc.

References basal, and IceBasalResistancePlasticLaw.drag().

Referenced by assembleSSAMatrix(), and correctBasalFrictionalHeating().

PetscScalar basalDragy	(	PetscScalar **	tauc,
		PISMVector2 **	uv,
		PetscInt	i,
		PetscInt	j
	)			const `[protected, virtual]`

Reimplemented in IceMISMIPModel.

Definition at line 36 of file iMbasal.cc.

References basal, and IceBasalResistancePlasticLaw.drag().

Referenced by assembleSSAMatrix(), and correctBasalFrictionalHeating().

PetscErrorCode basalSlidingHeatingSIA ( ) [protected, virtual]

Compute the basal sliding and frictional heating if (where) SIA sliding rule is used.

This routine is only called, by velocity(), if $\mu$ =mu_sliding is non-zero.

THIS KIND OF SIA SLIDING LAW IS A BAD IDEA. THAT'S WHY $\mu$ IS SET TO ZERO BY DEFAULT. See Appendix B of BBssasliding for the dangers in this mechanism.

This routine calls the SIA-type sliding law, which may return zero in the frozen base case; see basalVelocitySIA(). The basal sliding velocity is computed for all SIA points. This routine also computes the basal frictional heating. The basal velocity Vecs vub and vvb and the frictional heating Vec are fully over-written. Where the ice is floating, they all have value zero.

See correctBasalFrictionalHeating() for the SSA contribution.

Definition at line 409 of file iMsia.cc.

References basalVelocitySIA(), IceModelVec.begin_access(), config, IceGrid.dx, IceGrid.dy, EC, IceModelVec.end_access(), Enth3, NCConfigVariable.get(), IceModelVec2V.get_array(), IceModelVec2S.get_array(), EnthalpyConverter.getAbsTemp(), EnthalpyConverter.getPressureFromDepth(), IceModelVec3.getValZ(), grid, ice, IceModelVec2Mask.is_floating(), IceGrid.Lx, IceGrid.Ly, IceFlowLaw.rho, standard_gravity, PISMVector2.u, PISMVector2.v, vel_basal, vH, vMask, vRb, vWork2d, IceGrid.xm, IceGrid.xs, IceGrid.ym, and IceGrid.ys.

Referenced by velocity().

PetscScalar basalVelocitySIA	(	PetscScalar	x,
		PetscScalar	y,
		PetscScalar	H,
		PetscScalar	T,
		PetscScalar	alpha,
		PetscScalar	mu,
		PetscScalar	min_T
	)			const `[protected, virtual]`

Compute the coefficient of surface gradient, for basal sliding velocity as a function of driving stress in SIA regions.

THIS KIND OF SIA SLIDING LAW IS A BAD IDEA IN A THERMOMECHANICALLY-COUPLED MODEL. THAT'S WHY $\mu$ IS SET TO ZERO BY DEFAULT.

In SIA regions (= MASK_SHEET) a basal sliding law of the form

$\mathbf{U}_b = (u_b,v_b) = - C \nabla h$

is allowed. Here $\mathbf{U}_b$ is the horizontal velocity of the base of the ice (the "sliding velocity") and $h$ is the elevation of the ice surface. This procedure returns the positive coefficient $C$ in this relationship. This coefficient can depend of the thickness, the basal temperature, and the horizontal location.

The default version for IceModel here is location-independent pressure-melting-temperature-activated linear sliding. See Appendix B of BBssasliding for the dangers in this mechanism.

Parameter $\mu$ can be set by option -mu_sliding.

The returned coefficient is used in basalSlidingHeatingSIA().

This procedure is virtual and can be replaced by any derived class.

Reimplemented in IceEISModel, and IceCompModel.

Definition at line 682 of file iMsia.cc.

References IceFlowLaw.beta_CC_grad, ice, IceFlowLaw.rho, and standard_gravity.

Referenced by basalSlidingHeatingSIA().

PetscErrorCode bed_def_setup ( ) [protected, virtual]

Definition at line 21 of file iMbeddef.cc.

References beddef, IceGrid.com, config, NCConfigVariable.get_string(), grid, and PISMOptionsList().

Referenced by init_physics().

PetscErrorCode bed_def_step ( ) [protected, virtual]

Definition at line 64 of file iMbeddef.cc.

References beddef, config, dt, NCConfigVariable.get(), grid, last_bed_def_update, stdout_flags, PISMBedDef.update(), updateSurfaceElevationAndMask(), and IceGrid.year.

Referenced by step().

PetscErrorCode bootstrapFromFile ( const char * filename ) [virtual]

Read file and use heuristics to initialize PISM from typical 2d data available through remote sensing.

This procedure is called by the base class when option -boot_from is used.

See chapter 4 of the User's Manual. We read only 2D information from the bootstrap file.

Definition at line 34 of file iMbootstrap.cc.

References NCTool.close(), IceGrid.com, IceGrid.compute_vertical_levels(), config, NCTool.find_variable(), g, NCConfigVariable.get(), NCConfigVariable.get_flag(), PISMIO.get_grid_info(), grid, IceGrid.Lx, IceGrid.Ly, IceGrid.Lz, vnreport.nc, NCTool.open_for_reading(), PISMOptionsIsSet(), putTempAtDepth(), IceModelVec.range(), IceModelVec.regrid(), secpera, IceModelVec.set(), setMaskSurfaceElevation_bootstrap(), tau3, vbed, vbmr, verbPrintf(), vGhf, vH, vHmelt, vLatitude, vLongitude, vtillphi, and vuplift.

Referenced by set_vars_from_options().

PetscErrorCode bootstrapSetBedrockColumnTemp	(	PetscInt	i,
		PetscInt	j,
		PetscScalar	Ttopbedrock,
		PetscScalar	geothermflux,
		PetscScalar	bed_thermal_k
	)			`[protected, virtual]`

Set the temperatures in a column of bedrock based on a temperature at the top and a geothermal flux.

This procedure sets the temperatures in the bedrock that would be correct for our model in steady state. In steady state there would be a temperature at the top of the bed and a flux condition at the bottom and the temperatures would be linear in between.

Call Tb3.begin_access() before and Tb3.end_access() after this routine.

Definition at line 471 of file iMbootstrap.cc.

References grid, IceGrid.Mbz, IceModelVec3Bedrock.setInternalColumn(), Tb3, and IceGrid.zblevels.

Referenced by putTempAtDepth().

PetscErrorCode broadcastSSAVelocity ( bool updateVelocityAtDepth ) [protected, virtual]

At all SSA points, update the velocity field.

Once the vertically-averaged velocity field is computed by the SSA, this procedure updates the three-dimensional horizontal velocities $u$ and $v$ . Note that $w$ gets updated later by vertVelocityFromIncompressibility(). The three-dimensional velocity field is needed, for example, so that the temperature equation can include advection. Basal velocities also get updated.

Here is where the flag do_superpose controlled by option -super applies. If do_superpose is true then the just-computed velocity $v$ from the SSA is combined, in convex combination, to the stored velocity $u$ from the SIA computation:

$U = f(|v|)\, u + \left(1-f(|v|)\right)\, v.$

Here

$f(|v|) = 1 - (2/\pi) \arctan(10^{-4} |v|^2)$

is a function which decreases smoothly from 1 for $|v| = 0$ to 0 as $|v|$ becomes significantly larger than 100 m/a.

Definition at line 808 of file iMssa.cc.

References IceModelVec2Stag.begin_access(), IceModelVec.begin_access(), config, IceModelVec2Stag.end_access(), IceModelVec.end_access(), IceModelVec2V.get_array(), NCConfigVariable.get_flag(), IceModelVec3.getInternalColumn(), grid, MASK_DRAGGING_SHEET, MASK_SHEET, IceGrid.Mz, pi, secpera, PISMVector2.u, u3, uvbar, IceModelVec.v, PISMVector2.v, v3, IceModelVec2Mask.value(), vel_bar, vel_basal, vel_ssa, vMask, IceGrid.xm, IceGrid.xs, IceGrid.ym, and IceGrid.ys.

Referenced by IceExactSSAModel.run(), and velocity().

PetscErrorCode check_maximum_thickness ( ) [protected, virtual]

Check if the thickness of the ice is too large and extend the grid if necessary.

Extends the grid such that the new one has 2 (two) levels above the ice.

Definition at line 155 of file iMutil.cc.

References artm, IceModelVec.begin_access(), check_maximum_thickness_hook(), IceGrid.com, config, EC, IceModelVec.end_access(), Enth3, Enthnew3, IceModelVec3.extend_vertically(), NCConfigVariable.get_flag(), EnthalpyConverter.getEnthPermissive(), EnthalpyConverter.getPressureFromDepth(), grid, PISMSurfaceModel.ice_surface_temperature(), IceGrid.initial_Mz, Istag3, IceGrid.Lz, IceGrid.Mbz, IceGrid.Mz, vfnow.N, pism_filename_add_suffix(), IceModelVec.range(), save_snapshots, IceGrid.set_vertical_levels(), Sigma3, Sigmastag3, snapshots_file_is_ready, snapshots_filename, split_snapshots, surface, T3, tau3, taunew3, Tnew3, u3, v3, verbPrintf(), vH, vWork2d, w3, IceGrid.xm, IceGrid.xs, IceGrid.year, IceGrid.ym, IceGrid.ys, IceGrid.zblevels, and IceGrid.zlevels.

Referenced by massContExplicitStep().

PetscErrorCode check_maximum_thickness_hook ( const int old_Mz ) [protected, virtual]

Allows derived classes to extend their own IceModelVec3's in vertical.

Base class version does absolutely nothing.

Definition at line 282 of file iMutil.cc.

Referenced by check_maximum_thickness().

bool checkThinNeigh	(	PetscScalar	E,
		PetscScalar	NE,
		PetscScalar	N,
		PetscScalar	NW,
		PetscScalar	W,
		PetscScalar	SW,
		PetscScalar	S,
		PetscScalar	SE
	)			`[protected, virtual]`

Definition at line 654 of file iMtemp.cc.

Referenced by enthalpyAndDrainageStep(), and temperatureStep().

PetscErrorCode compute_bwp ( IceModelVec2S & result ) [protected, virtual]

Computes the subglacial (basal) water pressure.

$p_w = \alpha\, \frac{w}{w_{\text{max}}}\, \rho\, g\, H,$

where

$\alpha$ is the till pore water fraction (till_pw_fraction),
is the effective thickness of subglacial melt water (bwat)
$w_{\text{max}}$ is the maximum allowed value for (hmelt_max),
$\rho$ is the ice density (ice_density)
is the ice thickness (thk)

Result is set to invalid (_FillValue) where the ice is floating, there being no meaning to the above calculation.

Definition at line 831 of file iMreport.cc.

References IceModelVec.begin_access(), config, IceModelVec.end_access(), IceModelVec2Mask.fill_where_floating(), NCConfigVariable.get(), getBasalWaterPressure(), grid, pism_config_editor.result, IceModelVec.set_attr(), IceModelVec.set_attrs(), IceModelVec.set_name(), vbmr, vH, vHmelt, vMask, IceGrid.xm, IceGrid.xs, IceGrid.ym, and IceGrid.ys.

Referenced by compute_by_name().

PetscErrorCode compute_by_name	(	string	name,
		PetscScalar &	result
	)			`[protected, virtual]`

Compute a scalar diagnostic quantity by name.

Definition at line 1424 of file iMreport.cc.

References compute_ice_area(), compute_ice_area_floating(), compute_ice_area_grounded(), compute_ice_enthalpy(), compute_ice_volume(), config, dt, dvoldt, NCConfigVariable.get(), total_basal_ice_flux, total_sub_shelf_ice_flux, and total_surface_ice_flux.

PetscErrorCode compute_by_name	(	string	name,
		IceModelVec *&	result
	)			`[protected, virtual]`

Computes a diagnostic quantity given by name and returns a pointer to a pre-allocated work vector containing it.

For 2D quantities, result will point to vWork2d[0].

For 3D -- to Enthnew3, because we don't have a general-purpose 3D work vector.

Note that (depending on the quantity requested) vWork2d[1] might get used as a temporary storage.

Definition at line 1092 of file iMreport.cc.

References compute_bwp(), compute_cbar(), compute_cbase(), compute_cflx(), compute_csurf(), compute_cts(), compute_dhdt(), compute_enthalpybase(), compute_enthalpysurf(), compute_hardav(), compute_liqfrac(), compute_proc_ice_area(), compute_rank(), compute_taud(), compute_temp(), compute_temp_pa(), compute_tempbase(), compute_temppabase(), compute_tempsurf(), compute_uvelbase(), compute_uvelsurf(), compute_vvelbase(), compute_vvelsurf(), compute_wvelbase(), compute_wvelsurf(), Enthnew3, and vWork2d.

Referenced by update_viewers(), write_timeseries(), and write_variables().

PetscErrorCode compute_cbar ( IceModelVec2S & result ) [protected, virtual]

Computes cbar, the magnitude of vertically-integrated horizontal velocity of ice and masks out ice-free areas.

Definition at line 471 of file iMreport.cc.

References fill_missing.fill_value, IceModelVec2V.magnitude(), IceModelVec2S.mask_by(), secpera, IceModelVec.set_attr(), IceModelVec.set_attrs(), IceModelVec.set_glaciological_units(), IceModelVec.set_name(), vel_bar, vH, and IceModelVec.write_in_glaciological_units.

Referenced by compute_by_name().

PetscErrorCode compute_cbase	(	IceModelVec2S &	result,
		IceModelVec2S &	tmp
	)			`[protected, virtual]`

Computes cbase, the magnitude of horizontal velocity of ice at base of ice and masks out ice-free areas. Uses tmp as a preallocated temporary storage.

Definition at line 516 of file iMreport.cc.

References fill_missing.fill_value, IceModelVec3.getHorSlice(), IceModelVec2S.mask_by(), secpera, IceModelVec.set_attr(), IceModelVec.set_attrs(), IceModelVec.set_glaciological_units(), IceModelVec.set_name(), IceModelVec2S.set_to_magnitude(), u3, v3, vH, and IceModelVec.write_in_glaciological_units.

Referenced by compute_by_name().

PetscErrorCode compute_cflx	(	IceModelVec2S &	result,
		IceModelVec2S &	cbar
	)			`[protected, virtual]`

Computes cflx, the magnitude of vertically-integrated horizontal flux of ice.

Definition at line 493 of file iMreport.cc.

References fill_missing.fill_value, IceModelVec2S.mask_by(), IceModelVec.multiply_by(), secpera, IceModelVec.set_attr(), IceModelVec.set_attrs(), IceModelVec.set_glaciological_units(), IceModelVec.set_name(), vH, and IceModelVec.write_in_glaciological_units.

Referenced by compute_by_name().

PetscErrorCode compute_csurf	(	IceModelVec2S &	result,
		IceModelVec2S &	tmp
	)			`[protected, virtual]`

Computes csurf, the magnitude of horizontal velocity of ice at ice surface and masks out ice-free areas. Uses tmp as a preallocated temporary storage.

Definition at line 542 of file iMreport.cc.

References IceModelVec.begin_access(), IceModelVec.end_access(), fill_missing.fill_value, IceModelVec3.getSurfaceValues(), IceModelVec2S.mask_by(), secpera, IceModelVec.set_attr(), IceModelVec.set_attrs(), IceModelVec.set_glaciological_units(), IceModelVec.set_name(), IceModelVec2S.set_to_magnitude(), u3, v3, vH, and IceModelVec.write_in_glaciological_units.

Referenced by compute_by_name().

PetscErrorCode compute_cts ( IceModelVec3 & useForCTS ) [protected, virtual]

Compute the CTS field, CTS = E/E_s(p) and put in a global IceModelVec3 provided by user.

Definition at line 706 of file iMreport.cc.

References config, NCConfigVariable.get_flag(), IceModelVec.set(), IceModelVec.set_attrs(), IceModelVec.set_name(), and setCTSFromEnthalpy().

Referenced by compute_by_name().

PetscErrorCode compute_dhdt ( IceModelVec2S & result ) [protected, virtual]

Computes the rate of change of ice surface elevation as a sum of the bedrock uplift rate and the thickness rate of change.

Definition at line 870 of file iMreport.cc.

References IceModelVec.add(), IceModelVec.copy_from(), IceModelVec2S.mask_by(), IceModelVec.set_attrs(), IceModelVec.set_glaciological_units(), IceModelVec.set_name(), vdHdt, vH, vuplift, and IceModelVec.write_in_glaciological_units.

Referenced by compute_by_name().

PetscErrorCode compute_enthalpybase ( IceModelVec2S & result ) [protected, virtual]

Computes ice enthalpy at the base of ice.

Definition at line 951 of file iMreport.cc.

References Enth3, fill_missing.fill_value, IceModelVec3.getHorSlice(), IceModelVec2S.mask_by(), IceModelVec.set_attr(), IceModelVec.set_attrs(), IceModelVec.set_name(), and vH.

Referenced by compute_by_name().

PetscErrorCode compute_enthalpysurf ( IceModelVec2S & result ) [protected, virtual]

Computes ice enthalpy at the 1 m below the surface.

Definition at line 1049 of file iMreport.cc.

References IceModelVec.begin_access(), IceModelVec.end_access(), Enth3, fill_missing.fill_value, IceModelVec3.getSurfaceValues(), grid, pism_config_editor.result, IceModelVec.set_attr(), IceModelVec.set_attrs(), IceModelVec.set_name(), vH, IceGrid.xm, IceGrid.xs, IceGrid.ym, and IceGrid.ys.

Referenced by compute_by_name().

PetscErrorCode compute_hardav ( IceModelVec2S & result ) [protected, virtual]

Computes the vertically-averaged ice hardness.

Definition at line 779 of file iMreport.cc.

References IceFlowLaw.averagedHarness(), IceModelVec.begin_access(), compute_temp(), IceModelVec.end_access(), Enthnew3, IceFlowLaw.exponent(), IceModelVec3.getInternalColumn(), grid, ice, IceGrid.kBelowHeight(), pism_config_editor.result, IceModelVec.set_attr(), IceModelVec.set_attrs(), IceModelVec.set_name(), TEMPORARY_STRING_LENGTH, vH, IceGrid.xm, IceGrid.xs, IceGrid.ym, IceGrid.ys, and IceGrid.zlevels.

Referenced by compute_by_name().

PetscErrorCode compute_ice_area ( PetscScalar & result ) [protected, virtual]

Computes ice area, in m^2.

Definition at line 1286 of file iMreport.cc.

References IceModelVec.begin_access(), cell_area, IceGrid.com, IceGrid.dx, IceGrid.dy, IceModelVec.end_access(), grid, vH, IceModelVec.was_created(), IceGrid.xm, IceGrid.xs, IceGrid.ym, and IceGrid.ys.

Referenced by compute_by_name(), and massContExplicitStep().

PetscErrorCode compute_ice_area_floating ( PetscScalar & result ) [protected, virtual]

Computes floating ice area, in m^2.

Definition at line 1352 of file iMreport.cc.

References IceModelVec.begin_access(), cell_area, IceGrid.com, IceGrid.dx, IceGrid.dy, IceModelVec.end_access(), grid, IceModelVec2Mask.is_floating(), vH, vMask, IceModelVec.was_created(), IceGrid.xm, IceGrid.xs, IceGrid.ym, and IceGrid.ys.

Referenced by compute_by_name().

PetscErrorCode compute_ice_area_grounded ( PetscScalar & result ) [protected, virtual]

Computes grounded ice area, in m^2.

Definition at line 1318 of file iMreport.cc.

References IceModelVec.begin_access(), cell_area, IceGrid.com, IceGrid.dx, IceGrid.dy, IceModelVec.end_access(), grid, IceModelVec2Mask.is_grounded(), vH, vMask, IceModelVec.was_created(), IceGrid.xm, IceGrid.xs, IceGrid.ym, and IceGrid.ys.

Referenced by compute_by_name().

PetscErrorCode compute_ice_enthalpy ( PetscScalar & result ) [protected, virtual]

Computes the total ice enthalpy in J.

Units of the specific enthalpy field $E=$ (IceModelVec3.Enth3) are J kg-1. We integrate $E(t,x,y,z)$ over the entire ice fluid region $\Omega(t)$ , multiplying by the density to get units of energy:

$E_{\text{total}}(t) = \int_\Omega(t) E(t,x,y,z) \rho_i \,dx\,dy\,dz.$

Definition at line 1393 of file iMreport.cc.

References IceModelVec.begin_access(), IceGrid.com, config, IceGrid.dx, IceGrid.dy, IceModelVec.end_access(), Enth3, NCConfigVariable.get(), IceModelVec3.getInternalColumn(), grid, IceGrid.kBelowHeight(), vH, IceGrid.xm, IceGrid.xs, IceGrid.ym, IceGrid.ys, and IceGrid.zlevels.

Referenced by compute_by_name().

PetscErrorCode compute_ice_volume ( PetscScalar & result ) [protected, virtual]

Computes the ice volume, in m^3.

Definition at line 1254 of file iMreport.cc.

References IceModelVec.begin_access(), cell_area, IceGrid.com, IceGrid.dx, IceGrid.dy, IceModelVec.end_access(), grid, vH, IceModelVec.was_created(), IceGrid.xm, IceGrid.xs, IceGrid.ym, and IceGrid.ys.

Referenced by compute_by_name().

PetscErrorCode compute_liqfrac ( IceModelVec3 & useForLiqfrac ) [protected, virtual]

Compute the liquid fraction, and put in a global IceModelVec3 provided by user.

Definition at line 724 of file iMreport.cc.

References config, NCConfigVariable.get_flag(), IceModelVec.set(), IceModelVec.set_attrs(), IceModelVec.set_name(), and setLiquidFracFromEnthalpy().

Referenced by compute_by_name().

PetscErrorCode compute_proc_ice_area ( IceModelVec2S & result ) [protected, virtual]

Sets entrues of result to corresponding processor ranks.

Definition at line 678 of file iMreport.cc.

References IceModelVec.begin_access(), IceModelVec.end_access(), grid, pism_config_editor.result, IceModelVec.set_attrs(), IceModelVec.set_name(), IceModelVec.time_independent, vH, IceGrid.xm, IceGrid.xs, IceGrid.ym, and IceGrid.ys.

Referenced by compute_by_name().

PetscErrorCode compute_rank ( IceModelVec2S & result ) [protected, virtual]

Sets entrues of result to corresponding processor ranks.

Definition at line 662 of file iMreport.cc.

References IceModelVec.begin_access(), IceModelVec.end_access(), grid, IceGrid.rank, pism_config_editor.result, IceModelVec.set_attrs(), IceModelVec.set_name(), IceModelVec.time_independent, IceGrid.xm, IceGrid.xs, IceGrid.ym, and IceGrid.ys.

Referenced by compute_by_name().

PetscErrorCode compute_taud	(	IceModelVec2S &	result,
		IceModelVec2S &	tmp
	)			`[protected, virtual]`

Computes taud, the magnitude of driving shear stress at base of ice. Uses tmp as a preallocated temporary storage.

Definition at line 642 of file iMreport.cc.

References computeDrivingStress(), fill_missing.fill_value, IceModelVec2S.mask_by(), IceModelVec.set_attr(), IceModelVec.set_attrs(), IceModelVec.set_name(), IceModelVec2S.set_to_magnitude(), and vH.

Referenced by compute_by_name().

PetscErrorCode compute_temp ( IceModelVec3 & result ) [protected, virtual]

Compute the ice temperature corresponding to Enth3, and put in a given IceModelVec3.

Computes ice temperature from enthalpy (as a diagnostic quantity).

Definition at line 113 of file iMenthalpy.cc.

References IceModelVec.begin_access(), IceGrid.com, EC, IceModelVec.end_access(), Enth3, EnthalpyConverter.getAbsTemp(), IceModelVec3.getInternalColumn(), EnthalpyConverter.getPressureFromDepth(), grid, IceGrid.Mz, IceModelVec.set_attr(), IceModelVec.set_attrs(), IceModelVec.set_name(), vH, IceGrid.xm, IceGrid.xs, IceGrid.ym, IceGrid.ys, and IceGrid.zlevels.

Referenced by compute_by_name(), compute_hardav(), compute_tempbase(), and compute_tempsurf().

PetscErrorCode compute_temp_pa ( IceModelVec3 & result ) [protected, virtual]

Compute the pressure-adjusted temperature in degrees C corresponding to ice temperature, and put in a global IceModelVec3 provided by user.

Definition at line 742 of file iMreport.cc.

References IceModelVec.begin_access(), config, EC, IceModelVec.end_access(), Enth3, NCConfigVariable.get(), IceModelVec3.getInternalColumn(), EnthalpyConverter.getPATemp(), EnthalpyConverter.getPressureFromDepth(), grid, IceGrid.Mz, IceModelVec.set_attrs(), IceModelVec.set_name(), IceModelVec.shift(), vH, IceGrid.xm, IceGrid.xs, IceGrid.ym, IceGrid.ys, and IceGrid.zlevels.

Referenced by compute_by_name().

PetscErrorCode compute_tempbase ( IceModelVec2S & result ) [protected, virtual]

Computes ice temperature at the base of ice.

Definition at line 969 of file iMreport.cc.

References compute_temp(), Enthnew3, fill_missing.fill_value, IceModelVec3.getHorSlice(), IceModelVec2S.mask_by(), IceModelVec.set_attr(), IceModelVec.set_attrs(), IceModelVec.set_name(), and vH.

Referenced by compute_by_name().

PetscErrorCode compute_temppabase	(	IceModelVec3 &	hasPATemp,
		IceModelVec2S &	result
	)			`[protected, virtual]`

Computes pressure-adjusted ice temperature at the base of ice.

Definition at line 990 of file iMreport.cc.

References fill_missing.fill_value, IceModelVec3.getHorSlice(), IceModelVec2S.mask_by(), IceModelVec.set_attr(), IceModelVec.set_attrs(), IceModelVec.set_name(), and vH.

Referenced by compute_by_name().

PetscErrorCode compute_tempsurf ( IceModelVec2S & result ) [protected, virtual]

Computes ice temperature at the 1 m below the surface.

Definition at line 1011 of file iMreport.cc.

References IceModelVec.begin_access(), compute_temp(), IceModelVec.end_access(), Enthnew3, fill_missing.fill_value, IceModelVec3.getSurfaceValues(), grid, pism_config_editor.result, IceModelVec.set_attr(), IceModelVec.set_attrs(), IceModelVec.set_name(), vH, IceGrid.xm, IceGrid.xs, IceGrid.ym, and IceGrid.ys.

Referenced by compute_by_name().

PetscErrorCode compute_uvelbase ( IceModelVec2S & result ) [protected, virtual]

Definition at line 886 of file iMreport.cc.

References IceModelVec.begin_access(), IceModelVec.end_access(), fill_missing.fill_value, IceModelVec3.getHorSlice(), IceModelVec2S.mask_by(), secpera, IceModelVec.set_attr(), IceModelVec.set_attrs(), IceModelVec.set_glaciological_units(), IceModelVec.set_name(), u3, vH, and IceModelVec.write_in_glaciological_units.

Referenced by compute_by_name().

PetscErrorCode compute_uvelsurf ( IceModelVec2S & result ) [protected, virtual]

Computes uvelsurf, the x component of velocity of ice at ice surface.

Definition at line 570 of file iMreport.cc.

References IceModelVec.begin_access(), IceModelVec.end_access(), fill_missing.fill_value, IceModelVec3.getSurfaceValues(), IceModelVec2S.mask_by(), secpera, IceModelVec.set_attr(), IceModelVec.set_attrs(), IceModelVec.set_glaciological_units(), IceModelVec.set_name(), u3, vH, and IceModelVec.write_in_glaciological_units.

Referenced by compute_by_name().

PetscErrorCode compute_vvelbase ( IceModelVec2S & result ) [protected, virtual]

Definition at line 908 of file iMreport.cc.

References IceModelVec.begin_access(), IceModelVec.end_access(), fill_missing.fill_value, IceModelVec3.getHorSlice(), IceModelVec2S.mask_by(), secpera, IceModelVec.set_attr(), IceModelVec.set_attrs(), IceModelVec.set_glaciological_units(), IceModelVec.set_name(), v3, vH, and IceModelVec.write_in_glaciological_units.

Referenced by compute_by_name().

PetscErrorCode compute_vvelsurf ( IceModelVec2S & result ) [protected, virtual]

Computes vvelsurf, the y component of velocity of ice at ice surface.

Definition at line 593 of file iMreport.cc.

References IceModelVec.begin_access(), IceModelVec.end_access(), fill_missing.fill_value, IceModelVec3.getSurfaceValues(), IceModelVec2S.mask_by(), secpera, IceModelVec.set_attr(), IceModelVec.set_attrs(), IceModelVec.set_glaciological_units(), IceModelVec.set_name(), v3, vH, and IceModelVec.write_in_glaciological_units.

Referenced by compute_by_name().

PetscErrorCode compute_wvelbase ( IceModelVec2S & result ) [protected, virtual]

Definition at line 929 of file iMreport.cc.

References IceModelVec.begin_access(), IceModelVec.end_access(), fill_missing.fill_value, IceModelVec3.getHorSlice(), IceModelVec2S.mask_by(), secpera, IceModelVec.set_attr(), IceModelVec.set_attrs(), IceModelVec.set_glaciological_units(), IceModelVec.set_name(), vH, w3, and IceModelVec.write_in_glaciological_units.

Referenced by compute_by_name().

PetscErrorCode compute_wvelsurf ( IceModelVec2S & result ) [protected, virtual]

Computes wvelsurf, the vertical velocity of ice at ice surface.

Note that there is no need to mask out ice-free areas here, because wvelsurf is zero at those locations.

Definition at line 619 of file iMreport.cc.

References IceModelVec.begin_access(), IceModelVec.end_access(), fill_missing.fill_value, IceModelVec3.getSurfaceValues(), IceModelVec2S.mask_by(), secpera, IceModelVec.set_attr(), IceModelVec.set_attrs(), IceModelVec.set_glaciological_units(), IceModelVec.set_name(), vH, w3, and IceModelVec.write_in_glaciological_units.

Referenced by compute_by_name().

PetscErrorCode computeDrivingStress	(	IceModelVec2S &	vtaudx,
		IceModelVec2S &	vtaudy
	)			`[protected, virtual]`

Compute vector driving stress at base of ice on the regular grid.

Computes the driving stress at the base of the ice:

$\tau_d = - \rho g H \nabla h$

If use_eta is TRUE then the surface gradient $\nabla h$ is computed by the gradient of the transformed variable $\eta= H^{(2n+2)/n}$ (frequently, $\eta= H^{8/3}$ ). Because this quantity is more regular at ice sheet margins, we get a better surface gradient. When the thickness at a grid point is very small (below minThickEtaTransform in the procedure), the formula is slightly modified to give a lower driving stress.

In floating parts the surface gradient is always computed by the regular formula.

Saves it in user supplied Vecs vtaudx and vtaudy, which are treated as global. (I.e. we do not communicate ghosts.)

Definition at line 43 of file iMgeometry.cc.

References IceModelVec.begin_access(), config, IceModelVec2S.diff_x(), IceModelVec2S.diff_x_p(), IceModelVec2S.diff_y(), IceModelVec2S.diff_y_p(), IceGrid.dx, IceGrid.dy, IceModelVec.end_access(), IceFlowLaw.exponent(), NCConfigVariable.get_flag(), grid, ice, IceModelVec2Mask.is_grounded(), n, IceFlowLaw.rho, standard_gravity, vbed, vH, vh, vMask, IceGrid.xm, IceGrid.xs, IceGrid.ym, and IceGrid.ys.

Referenced by assembleSSARhs(), and compute_taud().

PetscErrorCode computeEffectiveViscosity	(	IceModelVec2S	vNuH[2],
		PetscReal	epsilon
	)			`[protected, virtual]`

Compute the product of the effective viscosity $\nu$ and ice thickness $H$ for the SSA model.

In PISM the product $\nu H$ can be

constant, or
can be computed with a constant ice hardness $\bar B$ (temperature-independent) but with dependence of the viscosity on the strain rates, or
it can depend on the strain rates and have a vertically-averaged ice hardness.

The flow law in ice stream and ice shelf regions must, for now, be a temperature-dependent Glen law. This is the only flow law we know how to convert to ``viscosity form''. (More general forms like Goldsby-Kohlstedt are not yet inverted.) The viscosity form is

$\nu(T^*,D) = \frac{1}{2} B(T^*) D^{(1/n)-1}\, D_{ij}$

where

$D_{ij} = \frac{1}{2} \left(\frac{\partial U_i}{\partial x_j} + \frac{\partial U_j}{\partial x_i}\right)$

is the strain rate tensor and $B$ is an ice hardness related to the ice softness $A(T^*)$ by

$B(T^*)=A(T^*)^{-1/n}$

in the case of a temperature dependent Glen-type law. (Here $T^*$ is the pressure-adjusted temperature.)

The effective viscosity is then

$\nu = \frac{\bar B}{2} \left[\left(\frac{\partial u}{\partial x}\right)^2 + \left(\frac{\partial v}{\partial y}\right)^2 + \frac{\partial u}{\partial x} \frac{\partial v}{\partial y} + \frac{1}{4} \left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}\right)^2 \right]^{(1-n)/(2n)}$

where in the temperature-dependent case

$\bar B = \frac{1}{H}\,\int_b^h B(T^*)\,dz$

This integral is approximately computed by the trapezoid rule.

In fact the integral is regularized as described in [SchoofStream]. The regularization constant $\epsilon$ is an argument to this procedure.

Also we put $\bar\nu H =$ constantNuHForSSA anywhere the ice is thinner than min_thickness_SSA. The geometry is not changed, but this has the effect of producing a shelf extension in ice free ocean, which affects the driving stress and the force balance at the calving front.

Definition at line 139 of file iMssa.cc.

References IceModelVec.begin_access(), IceModelVec.beginGhostComm(), IceGrid.com, config, IceGrid.dx, IceGrid.dy, EC, IceFlowLaw.effectiveViscosityColumn(), PolyThermalGPBLDIce.effectiveViscosityColumnFromEnth(), IceModelVec.end_access(), IceModelVec.endGhostComm(), Enth3, IceModelVec2V.get_array(), IceModelVec2S.get_array(), NCConfigVariable.get_flag(), EnthalpyConverter.getAbsTemp(), IceModelVec3.getInternalColumn(), EnthalpyConverter.getPressureFromDepth(), grid, ice, IceGrid.kBelowHeight(), leaveNuHAloneSSA, SSAStrengthExtension.min_thickness_for_extension(), IceGrid.Mz, SSAStrengthExtension.notional_strength(), IceModelVec.set(), ssaStrengthExtend, PISMVector2.u, PISMVector2.v, vel_ssa, vH, IceGrid.xm, IceGrid.xs, IceGrid.ym, IceGrid.ys, and IceGrid.zlevels.

Referenced by IceROSSModel.finishROSS(), and velocitySSA().

PetscErrorCode computeFlowUbarStats	(	PetscScalar *	gUbarmax,
		PetscScalar *	gUbarSIAav,
		PetscScalar *	gUbarstreamav,
		PetscScalar *	gUbarshelfav,
		PetscScalar *	gicegridfrac,
		PetscScalar *	gSIAgridfrac,
		PetscScalar *	gstreamgridfrac,
		PetscScalar *	gshelfgridfrac
	)			`[protected, virtual]`

Definition at line 27 of file iMreport.cc.

References MASK_DRAGGING_SHEET, and MASK_SHEET.

Referenced by summary().

PetscErrorCode computeMax2DSlidingSpeed ( ) [protected, virtual]

Compute the CFL constant associated to first-order upwinding for the sliding contribution to mass continuity.

This procedure computes the maximum horizontal speed in the SSA areas. In particular it computes CFL constant for the upwinding, in massContExplicitStep(), which applies to the basal component of mass flux.

That is, because the map-plane mass continuity is advective in the sliding case we have a CFL condition.

Definition at line 201 of file iMadaptive.cc.

References basal, IceModelVec.begin_access(), CFLmaxdt2D, IceGrid.com, config, IceGrid.dx, IceGrid.dy, IceModelVec.end_access(), NCConfigVariable.get(), IceModelVec2V.get_array(), NCConfigVariable.get_flag(), grid, IceModelVec2Mask.is_floating(), MASK_OCEAN_AT_TIME_0, secpera, IceModelVec2Mask.value(), vel_basal, vMask, IceGrid.xm, IceGrid.xs, IceGrid.ym, and IceGrid.ys.

Referenced by velocity().

PetscErrorCode computeMax3DVelocities ( ) [protected, virtual]

Compute the maximum velocities for time-stepping and reporting to user.

Computes the maximum magnitude of the components $u,v,w$ of the 3D velocity. Then sets CFLmaxdt, the maximum time step allowed under the Courant-Friedrichs-Lewy (CFL) condition on the horizontal advection scheme for age and for temperature.

Under BOMBPROOF there is no CFL condition for the vertical advection. The maximum vertical velocity is computed but it does not affect CFLmaxdt.

Definition at line 147 of file iMadaptive.cc.

References IceModelVec.begin_access(), CFLmaxdt, IceGrid.com, config, IceGrid.dx, IceGrid.dy, IceModelVec.end_access(), NCConfigVariable.get(), IceModelVec2S.get_array(), IceModelVec3.getInternalColumn(), gmaxu, gmaxv, gmaxw, grid, IceGrid.kBelowHeight(), secpera, u3, v3, vH, w3, IceGrid.xm, IceGrid.xs, IceGrid.ym, and IceGrid.ys.

Referenced by IceExactSSAModel.run(), and velocity().

PetscErrorCode computeMaxDiffusivity ( bool update_diffusivity_viewer ) [protected, virtual]

Compute the maximum diffusivity associated to the SIA deformational velocity.

The time-stepping scheme for mass continuity is explicit. It solves an equation which in the simple case of flat bed and the use of the nonsliding SIA is a pure diffusion,

$H_t = M - S - \nabla \cdot \mathbf{Q}$

where

$\mathbf{Q} = - D \nabla H$

is the horizontal ice flux and $D\ge 0$ is the diffusivity. Because the PDE is actually nonlinear, this diffusivity $D$ changes at every time step. The current procedure computes the maximum of the diffusivity on the grid.

Note that more generally the ice flow is driven by the driving stress $\tau_d = - \rho g H \nabla h$ which involves the surface slope, not the gradient of the thickness $H$ .

The time-stepping for the explicit scheme is controlled by equation (25) in [BBL], so that $\Delta t \sim \frac{\Delta x^2}{\max D}$ ; see also [MortonMayers]. But also

$\mathbf{Q} = \bar U\, H$

where $\bar U$ is the vertically-averaged horizontal velocity. Because of how the SIA calculation is currently factored, we compute $D$ here, for the purpose of adaptive time stepping, by the formula

$D = \frac{|\bar U| H}{|\nabla h|}.$

The potential division by zero in (harmless) areas of flat slope is avoided by addition of a constant to $\alpha = |\nabla h|$ .

See surfaceGradientSIA(), velocitySIAStaggered(), determineTimeStep(), and massContExplicitStep(); all of these are related calculations.

This method assumes IceModelVec2Stag uvbar holds correct deformational values of velocities coming out of the calculation in velocitySIAStaggered(). It also assumes the thickness in vH is up to date.

This method puts the maximum, over all staggered points, of the diffusivities into the global variable IceModel.gDmax.

If the user wants a run-time view of diffusivity, then that is updated here.

Definition at line 62 of file iMadaptive.cc.

References IceModelVec2Stag.begin_access(), IceModelVec.begin_access(), IceGrid.com, computeSIAVelocities, IceModelVec2Stag.end_access(), IceModelVec.end_access(), gDmax, IceModelVec2S.get_array(), grid, H0, IceModelVec.set(), IceModelVec.set_attrs(), IceModelVec.set_name(), surfaceGradientSIA(), uvbar, vH, IceModelVec2.view(), vWork2d, IceGrid.xm, IceGrid.xs, IceGrid.ym, and IceGrid.ys.

Referenced by determineTimeStep(), and summary().

PetscErrorCode computePhiFromBedElevation ( ) [protected, virtual]

Computes the till friction angle phi as a piecewise linear function of bed elevation, according to user options.

Computes the till friction angle $\phi(x,y)$ at a location, namely IceModel.vtillphi, as the following increasing, piecewise-linear function of the bed elevation $b(x,y)$ . Let

$M = (\phi_{\text{max}} - \phi_{\text{min}}) / (b_{\text{max}} - b_{\text{min}})$

be the slope of the nontrivial part. Then

$\phi(x,y) = \begin{cases} \phi_{\text{min}}, & b(x,y) \le b_{\text{min}}, \\ \phi_{\text{min}} + (b(x,y) - b_{\text{min}}) \,M, & b_{\text{min}} < b(x,y) < b_{\text{max}}, \\ \phi_{\text{max}}, & b_{\text{max}} \le b(x,y), \end{cases}$

The exception is if the point is marked as floating, in which case the till friction angle is set to the value phi_ocean.

The default values are vaguely suitable for Antarctica, perhaps:

phi_min = 5.0 degrees,
phi_max = 15.0 degrees,
topg_min = -1000.0 m,
topg_max = 1000.0 m,
phi_ocean = 10.0 degrees.

If the user gives option -topg_to_phi A,B,C,D then phi_ocean is not used. Instead, the same rule as above for grounded ice is used.

Definition at line 127 of file iMbasal.cc.

References IceModelVec.begin_access(), IceGrid.com, IceModelVec.end_access(), IceModelVec2S.get_array(), grid, IceModelVec2Mask.is_floating(), vbed, verbPrintf(), vMask, vtillphi, IceGrid.xm, IceGrid.xs, IceGrid.ym, and IceGrid.ys.

Referenced by initBasalTillModel().

PetscErrorCode correct_cell_areas ( ) [virtual]

Allocate and compute corrected cell areas. Uses linear interpolation to find latitudes and longitudes of grid corners, WGS84 parameters to compute cartesian coordinates of grid corners and vector products to compute areas of resulting triangles.

Note:: Latitude and longitude fields are not periodic, so computing corrected areas for cells at the grid boundary is not feasible. This should not matter, since these cells should be (and in most cases are) ice-free.; Using linear interpolation introduces errors in lon/lat coordinates of cell corners, but the corresponding ice volume error (present day Greenland, 5km grid) is about 3.11e-04 %.

Definition at line 434 of file iMutil.cc.

References PISMVars.add(), IceModelVec.begin_access(), cell_area, IceGrid.com, config, IceModelVec2S.create(), IceGrid.dx, IceGrid.dy, IceModelVec.end_access(), geo_x(), geo_y(), geo_z(), NCConfigVariable.get(), IceModelVec2S.get_array(), NCConfigVariable.get_flag(), grid, NCVariable.has(), mapping, IceGrid.Mx, IceGrid.My, IceModelVec.set(), IceModelVec.set_attrs(), IceModelVec.set_glaciological_units(), IceModelVec.time_independent, triangle_area(), variables, verbPrintf(), vLatitude, vLongitude, IceModelVec.write_in_glaciological_units, IceGrid.xm, IceGrid.xs, IceGrid.ym, and IceGrid.ys.

Referenced by misc_setup().

PetscErrorCode correctBasalFrictionalHeating ( ) [protected, virtual]

At SSA points, correct the previously-computed basal frictional heating.

Ice shelves have zero basal friction heating.

Definition at line 884 of file iMssa.cc.

References basalDragx(), basalDragy(), IceModelVec.begin_access(), config, IceModelVec.end_access(), IceModelVec2S.get_array(), IceModelVec2V.get_array(), NCConfigVariable.get_flag(), grid, IceModelVec2Mask.is_floating(), MASK_DRAGGING_SHEET, PISMVector2.u, PISMVector2.v, IceModelVec2Mask.value(), vel_basal, vMask, vRb, vtauc, IceGrid.xm, IceGrid.xs, IceGrid.ym, and IceGrid.ys.

Referenced by velocity().

PetscErrorCode correctSigma ( ) [protected, virtual]

At SSA points, correct the previously-computed volume strain heating (dissipation heating).

Documented in BBssasliding.

Definition at line 926 of file iMssa.cc.

References IceModelVec.begin_access(), IceFlowLaw.beta_CC_grad, config, IceGrid.dx, IceGrid.dy, e, EC, IceModelVec.end_access(), Enth3, IceFlowLaw.exponent(), NCConfigVariable.get(), IceModelVec2V.get_array(), IceModelVec2S.get_array(), NCConfigVariable.get_flag(), EnthalpyConverter.getAbsTemp(), IceModelVec3.getInternalColumn(), EnthalpyConverter.getPressureFromDepth(), grid, IceFlowLaw.hardnessParameter(), ice, IceGrid.kBelowHeight(), MASK_DRAGGING_SHEET, MASK_SHEET, IceGrid.Mz, pi, secpera, Sigma3, PISMVector2.u, PISMVector2.v, IceModelVec2Mask.value(), vel_ssa, vH, vMask, IceGrid.xm, IceGrid.xs, IceGrid.ym, IceGrid.ys, and IceGrid.zlevels.

Referenced by velocity().

PetscErrorCode countCFLViolations ( PetscScalar * CFLviol ) [protected, virtual]

Because of the -skip mechanism it is still possible that we can have CFL violations: count them.

This applies to the horizontal part of the three-dimensional advection problem solved by IceModel.ageStep() and the advection, ice-only part of the problem solved by temperatureStep(). These methods use a fine vertical grid, and so we consider CFL violations on that same fine grid. (FIXME: should we actually use the fine grid?)

Communication is needed to determine total CFL violation count over entire grid. It is handled by temperatureAgeStep(), not here.

Definition at line 345 of file iMadaptive.cc.

References IceModelVec.begin_access(), dtTempAge, IceGrid.dx, IceGrid.dy, IceModelVec.end_access(), IceModelVec3.getInternalColumn(), grid, IceGrid.kBelowHeight(), u3, v3, vH, IceGrid.xm, IceGrid.xs, IceGrid.ym, and IceGrid.ys.

Referenced by energyStep().

PetscErrorCode create_timeseries ( ) [protected]

Creates DiagnosticTimeseries objects used to store and report scalar diagnostic quantities.

Definition at line 113 of file iMtimeseries.cc.

References grid, DiagnosticTimeseries.output_filename, Timeseries.set_attr(), Timeseries.set_dimension_units(), Timeseries.set_units(), timeseries, ts_filename, and ts_vars.

Referenced by init_timeseries().

PetscErrorCode createVecs ( ) [protected, virtual]

Allocate all IceModelVecs defined in IceModel.

This procedure allocates the memory used to store model state, diagnostic and work vectors and sets metadata.

Default values should not be set here; please use set_vars_from_options().

All the memory allocated here is freed by IceModelVecs' destructors.

Reimplemented in IceEISModel, IceROSSModel, IceMISMIPModel, IceUnitModel, IceCompModel, and IceExactSSAModel.

Definition at line 116 of file iceModel.cc.

References acab, PISMVars.add(), artm, IceGrid.com, config, IceModelVec2Stag.create(), IceModelVec2V.create(), IceModelVec2Mask.create(), IceModelVec2S.create(), IceModelVec3Bedrock.create(), IceModelVec3.create(), Enth3, NCConfigVariable.get(), NCConfigVariable.get_flag(), grid, MASK_DRAGGING_SHEET, MASK_FLOATING, MASK_ICE_FREE_BEDROCK, MASK_ICE_FREE_OCEAN, MASK_OCEAN_AT_TIME_0, MASK_SHEET, MY_WIDE_STENCIL, IceModelVec.output_data_type, secpera, IceModelVec.set_attr(), IceModelVec.set_attrs(), IceModelVec.set_glaciological_units(), shelfbmassflux, shelfbtemp, Sigma3, T3, tau3, Tb3, IceModelVec.time_independent, u3, uvbar, v3, variables, vbed, vbmr, vdHdt, vel_bar, vel_basal, vel_ssa, verbPrintf(), vGhf, vH, vh, vHmelt, vLatitude, vLongitude, vMask, vRb, vtauc, vtillphi, vuplift, w3, and IceModelVec.write_in_glaciological_units.

Referenced by init().

PetscErrorCode deallocate_internal_objects ( ) [protected, virtual]

De-allocate internal objects.

This includes Vecs that are not in an IceModelVec, SSA tools and the bed deformation model.

Definition at line 417 of file iceModel.cc.

References destroySSAobjects().

Referenced by ~IceModel().

PetscErrorCode destroySSAobjects ( ) [protected, virtual]

Deallocate SSA tools.

Definition at line 60 of file iMssa.cc.

References SSADA, SSAKSP, SSARHS, SSAStiffnessMatrix, and SSAX.

Referenced by deallocate_internal_objects().

PetscErrorCode determineTimeStep ( const bool doTemperatureCFL ) [protected, virtual]

Use various stability criteria to determine the time step for an evolution run.

The main loop in run() approximates many physical processes. Several of these approximations, including the mass continuity and temperature equations in particular, involve stability criteria. This procedure builds the length of the next time step by using these criteria and by incorporating choices made by options (e.g. -max_dt) and by derived classes.

Definition at line 276 of file iMadaptive.cc.

References adaptReasonFlag, adaptTimeStepDiffusivity(), CFLmaxdt, CFLmaxdt2D, computeMaxDiffusivity(), computeSIAVelocities, config, doAdaptTimeStep, dt, dt_force, dt_from_cfl, IceGrid.end_year, NCConfigVariable.get(), NCConfigVariable.get_flag(), grid, maxdt_temporary, secpera, skipCountDown, view_diffusivity, and IceGrid.year.

Referenced by step().

PetscErrorCode diffuseHmelt ( ) [protected, virtual]

Apply explicit time step for pure diffusion to basal layer of melt water.

See preprint BBssasliding .

Uses vWork2d[0] to temporarily store new values for Hmelt.

Definition at line 386 of file iMbasal.cc.

References IceModelVec.beginGhostComm(), config, dtTempAge, IceGrid.dx, IceGrid.dy, IceModelVec.end_access(), IceModelVec.endGhostComm(), NCConfigVariable.get(), IceModelVec2S.get_array(), grid, L, secpera, vHmelt, vWork2d, IceGrid.xm, IceGrid.xs, IceGrid.ym, and IceGrid.ys.

Referenced by step().

PetscErrorCode dumpToFile ( const char * filename ) [protected, virtual]

Definition at line 68 of file iMIO.cc.

References NCTool.append_time(), NCTool.close(), config, dt, global_attributes, grid, mapping, vnreport.nc, ocean, PISMIO.open_for_writing(), overrides, PISMOptionsIsSet(), secpera, surface, NCConfigVariable.update_from(), NCGlobalAttributes.write(), NCConfigVariable.write(), write_extra_fields(), PISMComponent.write_model_state(), PISMSurfaceModel.write_model_state(), write_model_state(), and IceGrid.year.

Referenced by endOfTimeStepHook(), and writeFiles().

int endOfTimeStepHook ( ) [protected, virtual]

Catch signals -USR1 and -TERM; in the former case save and continue; in the latter, save and stop.

Signal SIGTERM makes PISM end, saving state under original -o name (or default name). We also add an indication to the history attribute of the output NetCDF file.

Signal SIGUSR1 makes PISM save state under a filename based on the the name of the executable (e.g. pismr or pismv) and the current model year. In addition the time series (-ts_file, etc.) is flushed out There is no indication of these actions in the history attribute of the output (-o) NetCDF file because there is no effect on it, but there is an indication at stdout.

Signal SIGUSR2 makes PISM flush time-series, without saving model state.

Definition at line 50 of file iMutil.cc.

References IceGrid.com, dumpToFile(), executable_short_name, grid, pism_signal, stampHistory(), TEMPORARY_STRING_LENGTH, timeseries, verbPrintf(), and IceGrid.year.

Referenced by run().

PetscErrorCode energyStats	(	PetscScalar	iarea,
		bool	useHomoTemp,
		PetscScalar &	gmeltfrac,
		PetscScalar &	gtemp0
	)			`[protected, virtual]`

Computes fraction of the base which is melted and the ice basal temperature at the center of the ice sheet.

Communication occurs here.

Definition at line 157 of file iMreport.cc.

References IceModelVec.begin_access(), IceGrid.com, IceGrid.dx, IceGrid.dy, e, EC, IceModelVec.end_access(), Enth3, IceModelVec2S.get_array(), EnthalpyConverter.getAbsTemp(), IceModelVec3.getHorSlice(), EnthalpyConverter.getPressureFromDepth(), grid, EnthalpyConverter.isTemperate(), IceGrid.Mx, IceGrid.My, vH, vWork2d, IceGrid.xm, IceGrid.xs, IceGrid.ym, and IceGrid.ys.

Referenced by summary().

PetscErrorCode energyStep ( ) [protected, virtual]

Manage the solution of the energy equation, and related parallel communication.

Calls the method enthalpyAndDrainageStep(), or temperatureStep() if do_cold_ice_methods == true.

This method (energyStep()) must update these four fields:

IceModelVec3 Enth3
IceModelVec3Bedrock Tb3
IceModelVec2 vbasalMeltRate
IceModelVec2 vHmelt That is, energyStep() is in charge of calling other methods that update, and then it is in charge of doing the ghost communication as needed. If do_cold_ice_methods == true, then energyStep() must also update this field
IceModelVec3 T3

Definition at line 44 of file iMtemp.cc.

References IceModelVec3.beginGhostCommTransfer(), CFLviolcount, IceGrid.com, config, countCFLViolations(), IceModelVec3.endGhostCommTransfer(), Enth3, enthalpyAndDrainageStep(), Enthnew3, NCConfigVariable.get_flag(), getVerbosityLevel(), grid, IceGrid.Mbz, IceGrid.Mx, IceGrid.My, IceGrid.Mz, setEnth3FromT3_ColdIce(), stdout_flags, T3, temperatureStep(), Tnew3, and verbPrintf().

Referenced by step().

PetscErrorCode enthalpyAndDrainageStep	(	PetscScalar *	vertSacrCount,
		PetscScalar *	liquifiedVol
	)			`[protected, virtual]`

Update enthalpy field based on conservation of energy in ice and bedrock.

This method is documented by the page BOMBPROOF, PISM's numerical scheme for conservation of energy.

This method uses instances of combinedSystemCtx, bedrockOnlySystemCtx, and iceenthOnlySystemCtx.

This method modifies IceModelVec3 Enthnew3, IceModelVec3Bedrock Tb3, IceModelVec2S vBasalMeltRate, and IceModelVec2S vHmelt. No communication of ghosts is done for any of these fields.

Regarding drainage, we move the liquid water fraction which is in excess of the fixed fraction omega_max in a column segment [z,z+dz] to the base. Heuristic: Once liquid water fraction exceeds a cap, all of it goes to the base. Follows [Greve97Greenland] and references therein.

Definition at line 336 of file iMenthalpy.cc.

References artm, IceModelVec.begin_access(), IceFlowLaw.c_p, checkThinNeigh(), IceGrid.com, config, copyColumn(), double, dtTempAge, IceGrid.dx, IceGrid.dy, IceGrid.dz_fine, EC, IceModelVec.end_access(), combinedSystemCtx.Enth, iceenthOnlySystemCtx.Enth, Enth3, combinedSystemCtx.Enth_s, iceenthOnlySystemCtx.Enth_s, Enthnew3, bedrockOnlySystemCtx.extractHeatFluxFromSoln(), NCConfigVariable.get(), NCConfigVariable.get_flag(), EnthalpyConverter.getAbsTemp(), getEnthalpyCTSColumn(), EnthalpyConverter.getEnthAtWaterFraction(), EnthalpyConverter.getEnthPermissive(), EnthalpyConverter.getMeltingTemp(), IceModelVec3.getPlaneStar(), EnthalpyConverter.getPressureFromDepth(), IceModelVec3.getValColumn(), IceModelVec3Bedrock.getValColumnPL(), getVerbosityLevel(), EnthalpyConverter.getWaterFractionLimited(), grid, ice, PISMSurfaceModel.ice_surface_temperature(), planeStar.ij, planeStar.im1, iceenthOnlySystemCtx.initAllColumns(), bedrockOnlySystemCtx.initAllColumns(), combinedSystemCtx.initAllColumns(), planeStar.ip1, IceModelVec2Mask.is_floating(), EnthalpyConverter.isLiquified(), issounding(), planeStar.jm1, planeStar.jp1, IceFlowLaw.k, L, IceGrid.Mbz_fine, IceGrid.Mz_fine, ocean, check_stationarity.p, PISMOptionsIsSet(), reportColumn(), reportColumnSolveError(), IceFlowLaw.rho, secpera, iceenthOnlySystemCtx.setBoundaryValuesThisColumn(), combinedSystemCtx.setBoundaryValuesThisColumn(), bedrockOnlySystemCtx.setBoundaryValuesThisColumn(), IceModelVec3Bedrock.setColumn(), IceModelVec3.setColumn(), columnSystemCtx.setIndicesAndClearThisColumn(), iceenthOnlySystemCtx.setLevel0EqnThisColumn(), iceenthOnlySystemCtx.setSchemeParamsThisColumn(), combinedSystemCtx.setSchemeParamsThisColumn(), IceModelVec3.setValColumnPL(), IceModelVec3Bedrock.setValColumnPL(), PISMOceanModel.shelf_base_mass_flux(), PISMOceanModel.shelf_base_temperature(), shelfbmassflux, shelfbtemp, iceenthOnlySystemCtx.Sigma, combinedSystemCtx.Sigma, Sigma3, iceenthOnlySystemCtx.solveThisColumn(), combinedSystemCtx.solveThisColumn(), bedrockOnlySystemCtx.solveThisColumn(), surface, combinedSystemCtx.Tb, bedrockOnlySystemCtx.Tb, Tb3, iceenthOnlySystemCtx.u, combinedSystemCtx.u, u3, updateHmelt, iceenthOnlySystemCtx.v, combinedSystemCtx.v, v3, vbmr, vGhf, vH, vHmelt, iceenthOnlySystemCtx.viewConstants(), bedrockOnlySystemCtx.viewConstants(), combinedSystemCtx.viewConstants(), EnthalpyConverter.viewConstants(), vMask, vRb, combinedSystemCtx.w, iceenthOnlySystemCtx.w, w3, IceGrid.xm, IceGrid.xs, IceGrid.year, IceGrid.ym, IceGrid.ys, and IceGrid.zlevels_fine.

Referenced by energyStep().

PetscErrorCode excessToFromBasalMeltLayer	(	const PetscScalar	rho_c,
		const PetscScalar	z,
		const PetscScalar	dz,
		PetscScalar *	Texcess,
		PetscScalar *	Hmelt
	)			`[protected, virtual]`

Compute the melt water which should go to the base if $T$ is above pressure-melting.

Definition at line 493 of file iMtemp.cc.

References allowAboveMelting, IceGrid.dx, IceGrid.dy, grid, ice, IceFlowLaw.latentHeat, IceFlowLaw.rho, and updateHmelt.

Referenced by temperatureStep().

PetscErrorCode extras_max_timestep	(	double	t_years,
		double &	dt_years
	)			`[protected]`

Computes the maximum time-step we can take and still hit all the requested years.

Sets dt_years to -1 if any time-step is OK.

Definition at line 512 of file iMtimeseries.cc.

References config, extra_times, NCConfigVariable.get_flag(), and save_extra.

Referenced by step().

PetscScalar getBasalWaterPressure	(	PetscScalar	thk,
		PetscScalar	bwat,
		PetscScalar	bmr,
		PetscScalar	frac,
		PetscScalar	hmelt_max
	)			const `[protected, virtual]`

Compute modeled pressure in subglacial liquid water using thickness of water layer, and possibly the melt rate, at the base.

This procedure provides a simple model of basal water pressure $p_w$ , which is modeled as a function of the thickness of the basal stored water plus (optionally) the basal melt rate.

Input bwat is thickness of basal water. Input bmr is the basal melt rate. Because both bwat and bmr are zero at points where base of ice is below the pressure-melting temperature, the modeled basal water pressure is zero when the base is frozen.

The inequality bwat $\le$ hmelt_max is required at input, and an error is thrown if not.

Regarding the physics, compare the water pressure computed by formula (4) in [Ritzetal2001], where the pressure is a function of sea level and bed elevation. Also, the method using "elevation of the bed at the grounding line" as in [LingleBrown1987] is not implementable because that elevation is at an unknowable location. (We are not doing a flow line model!)

Several options control the water pressure model:

-[no_]bmr_enhance toggle the basal melt rate dependency in water pressure; DEFAULT IS OFF
-bmr_enhance_scale sets the value; defaults to 0.10, since 10 cm/a is a significant enough level of basal melt rate to cause weakening effect for saturated till; argument is in m/a; must set -bmr_enhance for this to have effect
-plastic_pwfrac controls parameter till_pw_fraction
-[no_]thk_eff toggle the thickness effect: for smaller thicknesses there is a reduction in basal water pressure, a conceptual near-margin drainage mechanism; DEFAULT IS OFF
-thk_eff_reduced factor by which basal water pressure is reduced; default is 0.97; must set -thk_eff for this to have effect
-thk_eff_H_high maximum thickness at which effect is applied; default is 2000 m; must set -thk_eff for this to have any effect
-thk_eff_H_low thickness at which thickness effect is full strength; default is 1000 m; must set -thk_eff for this to have any effect

At several places in the code the effective pressure on the mineral part of the till is computed by these lines, which are recommended for this purpose:

p_over = ice->rho * standard_gravity * thk; // the pressure of the weight of the ice

p_eff = p_over - getBasalWaterPressure(thk, bwat, bmr, frac, hmelt_max);

Definition at line 247 of file iMbasal.cc.

References IceGrid.com, config, NCConfigVariable.get(), NCConfigVariable.get_flag(), grid, ice, IceFlowLaw.rho, and standard_gravity.

Referenced by compute_bwp(), and updateYieldStressUsingBasalWater().

PetscErrorCode getEnthalpyCTSColumn	(	PetscScalar	p_air,
		PetscScalar	thk,
		PetscInt	ks,
		const PetscScalar *	Enth,
		const PetscScalar *	w,
		PetscScalar *	lambda,
		PetscScalar **	Enth_s
	)			`[protected, virtual]`

Compute the CTS value of enthalpy in an ice column, and the lambda for BOMBPROOF.

Return argument Enth_s[Mz] has the enthalpy value for the pressure-melting temperature at the corresponding z level.

Parameters:

	p_air	atmospheric pressure
	thk	ice thickness
	ks	index of the level just below the surface

Definition at line 224 of file iMenthalpy.cc.

References IceFlowLaw.c_p, IceGrid.dz_fine, EC, EnthalpyConverter.getEnthalpyCTS(), EnthalpyConverter.getPressureFromDepth(), grid, ice, IceFlowLaw.k, IceGrid.Mz_fine, IceFlowLaw.rho, secpera, and IceGrid.zlevels_fine.

Referenced by enthalpyAndDrainageStep().

virtual IceFlowLaw* getIceFlowLaw ( ) [virtual]

Definition at line 84 of file iceModel.hh.

References ice.

Referenced by main().

virtual IceFlowLawFactory& getIceFlowLawFactory ( ) [virtual]

Definition at line 83 of file iceModel.hh.

References iceFactory.

PetscScalar grainSizeVostok ( PetscScalar age ) const [protected, virtual]

Use the Vostok core as a source of a relationship between the age of the ice and the grain size.

A data set is interpolated here. The intention is that the softness of the ice has nontrivial dependence on its age, through its grainsize, because of variable dustiness of the global climate. The grainsize is partly determined by at which point in the glacial cycle the given ice fell as snow.

The data is from DeLaChapelleEtAl98 and LipenkovEtAl89 . In particular, Figure A2 in the former reference was hand-sampled with an attempt to include the ``wiggles'' in that figure. Ages of the oldest ice (>= 300 ka) were estimated in a necessarily ad hoc way. The age value of 10000 ka was added simply to give interpolation for very old ice; ages beyond that get constant extrapolation. Linear interpolation is done between the samples.

Definition at line 38 of file iMgrainsize.cc.

References IceGrid.com, grid, and secpera.

Referenced by velocitySIAStaggered().

PetscErrorCode grid_setup ( ) [virtual]

Sets up the computational grid.

There are two cases here:

1) Initializing from a PISM ouput file, in which case all the options influencing the grid (currently: -Mx, -My, -Mz, -Mbz, -Lx, -Ly, -Lz, -z_spacing, -zb_spacing) are ignored.

2) Initializing using defaults, command-line options and (possibly) a bootstrapping file. Derived classes requiring special grid setup should reimplement IceGrid.set_grid_from_options().

No memory allocation should happen here.

Definition at line 225 of file iMinit.cc.

References NCTool.close(), IceGrid.com, IceGrid.createDA(), pism_config_editor.filename, NCTool.find_variable(), NCTool.get_att_text(), PISMIO.get_grid(), grid, ignore_option(), mapping, vnreport.nc, NCTool.open_for_reading(), PISMOptionsString(), NCConfigVariable.print(), NCConfigVariable.read(), set_grid_defaults(), set_grid_from_options(), set_time_from_options(), IceGrid.start_year, verbPrintf(), and IceGrid.year.

PetscErrorCode horizontalVelocitySIARegular ( ) [protected, virtual]

Update regular grid horizontal velocities u,v at depth for SIA regions.

The procedure velocitySIAStaggered() computes several scalar quantities at depth (the details of which are too complicated to explain). These quantities correspond to three-dimensional arrays. This procedure takes those quantities and computes the three-dimensional arrays for the horizontal components $u$ and $v$ of the velocity field.

The vertical component $w$ of the velocity field is computed later by vertVelocityFromIncompressibility().

Definition at line 596 of file iMsia.cc.

References IceModelVec.begin_access(), config, IceModelVec.end_access(), NCConfigVariable.get(), IceModelVec2S.get_array(), IceModelVec3.getInternalColumn(), grid, Istag3, IceGrid.Mz, IceModelVec3.setInternalColumn(), u3, IceModelVec.v, v3, vel_basal, vWork2d, IceGrid.xm, IceGrid.xs, IceGrid.ym, and IceGrid.ys.

Referenced by velocity().

PetscErrorCode ice_mass_bookkeeping ( ) [protected, virtual]

Computes total ice fluxes in kg s-1 at 3 interfaces:

the ice-atmosphere interface: gets surface mass balance rate from PISMSurfaceModel *surface,
the ice-ocean interface at the bottom of ice shelves: gets ocean-imposed basal melt rate from PISMOceanModel *ocean, and
the ice-bedrock interface: gets basal melt rate from IceModelVec2S vbmr.

A unit-conversion occurs for all three quantities, from ice-equivalent m s-1 to kg s-1. The sign convention about these fluxes is that positve flux means ice is being added to the ice fluid volume at that interface.

These quantities should be understood as instantaneous at the beginning of the time-step. Multiplying by dt will not necessarily give the corresponding change from the beginning to the end of the time-step.

FIXME: The calving rate can be computed by post-processing: divoldt = surface_ice_flux * iarea + basal_ice_flux * iareag + sub_shelf_ice_flux * iareaf + calving_flux_vol_rate

Definition at line 1513 of file iMreport.cc.

References acab, IceModelVec.begin_access(), cell_area, IceGrid.com, config, dt, IceGrid.dx, IceGrid.dy, IceModelVec.end_access(), NCConfigVariable.get(), NCConfigVariable.get_flag(), grid, PISMSurfaceModel.ice_surface_mass_flux(), IceModelVec2Mask.is_grounded(), MASK_FLOATING, ocean, secpera, PISMOceanModel.shelf_base_mass_flux(), shelfbmassflux, surface, total_basal_ice_flux, total_sub_shelf_ice_flux, total_surface_ice_flux, IceModelVec2Mask.value(), vbmr, vH, vMask, IceModelVec.was_created(), IceGrid.xm, IceGrid.xs, IceGrid.year, IceGrid.ym, and IceGrid.ys.

Referenced by step().

PetscErrorCode init ( )

Manage the initialization of the IceModel object.

Please see the documenting comments of the functions called below to find explanations of their intended uses.

The IceModel initialization sequence is this:

1) Initialize the computational grid:

2) Process the options:

3) Memory allocation:

4) Initialize the IceFlowLaw and (possibly) other physics.

5) Initialize atmosphere and ocean couplers:

6) Fill the model state variables (from a PISM output file, from a bootstrapping file using some modeling choices or using formulas). Calls IceModel.regrid()

7) Report grid parameters:

8) Allocate SSA tools and work vectors:

9) Miscellaneous stuff: set up the bed deformation model, initialize the basal till model, initialize snapshots. This has to happen *after* regridding.

The following flow-chart illustrates the process. initialization-sequence.dot

IceModel initialization sequence

Definition at line 714 of file iceModel.cc.

References allocate_internal_objects(), IceGrid.com, createVecs(), grid, init_couplers(), init_physics(), misc_setup(), model_state_setup(), pism_wait_for_gdb(), PISMOptionsIsSet(), report_grid_parameters(), and setFromOptions().

Referenced by main().

PetscErrorCode init_couplers ( ) [virtual]

Initializes atmosphere and ocean couplers.

Reimplemented in IceEISModel, and IceMISMIPModel.

Definition at line 464 of file iMinit.cc.

References IceGrid.com, grid, PISMComponent.init(), PISMSurfaceModel.init(), ocean, surface, variables, and verbPrintf().

Referenced by init().

PetscErrorCode init_extras ( ) [protected]

Initialize the code saving spatially-variable diagnostic quantities.

Definition at line 314 of file iMtimeseries.cc.

References IceGrid.com, config, current_extra, ends_with(), extra_file_is_ready, extra_filename, extra_times, extra_vars, NCConfigVariable.get_flag(), PISMVars.get_variables(), grid, parse_times(), PISMOptionsIsSet(), PISMOptionsString(), save_extra, split_extra, variables, and verbPrintf().

Referenced by misc_setup().

PetscErrorCode init_physics ( ) [virtual]

Initialize some physical parameters.

This is the place for all non-trivial initialization of physical parameters. ("Non-trivial" means that the initialization requires more than just setting a value of a parameter. Such trivial changes can go here too, or earlier.) Also, this is the good place to set those parameters that a user should not be able to override using a command-line option.

This method is called after memory allocation but before filling any of IceModelVecs because all the physical parameters should be initialized before setting up the coupling or filling model-state variables.

In the base class IceModel we just initialize the IceFlowLaw and the EnthalpyConverter.

Reimplemented in IceEISModel, IcePSTexModel, IceROSSModel, IceMISMIPModel, IceCompModel, and IceExactSSAModel.

Definition at line 407 of file iMinit.cc.

References bed_def_setup(), IceGrid.com, config, IceFlowLawFactory.create(), EC, NCConfigVariable.get_flag(), grid, ice, ICE_GPBLD, iceFactory, IceFlowLaw.printInfo(), IceFlowLaw.setFromOptions(), IceFlowLawFactory.setFromOptions(), IceFlowLawFactory.setType(), and verbPrintf().

Referenced by init().

PetscErrorCode init_snapshots ( ) [protected]

Initializes the snapshot-saving mechanism.

Definition at line 416 of file iMIO.cc.

References IceGrid.com, current_snapshot, ends_with(), grid, parse_times(), PISMOptionsIsSet(), PISMOptionsString(), save_snapshots, set_output_size(), snapshot_times, snapshot_vars, snapshots_file_is_ready, snapshots_filename, split_snapshots, and verbPrintf().

Referenced by misc_setup().

PetscErrorCode init_timeseries ( ) [protected]

Initializes the code writing scalar time-series.

Definition at line 25 of file iMtimeseries.cc.

References NCTool.close(), IceGrid.com, config, create_timeseries(), current_ts, NCConfigVariable.get_string(), grid, vnreport.nc, PISMIO.open_for_writing(), parse_times(), PISMOptionsIsSet(), PISMOptionsString(), save_ts, ts_filename, ts_times, ts_vars, and verbPrintf().

Referenced by misc_setup().

PetscErrorCode init_viewers ( ) [protected, virtual]

Initialize run-time diagnostic viewers.

Definition at line 144 of file iMviewers.cc.

References IceGrid.com, config, NCConfigVariable.get(), grid, IceGrid.Lz, map_viewers, NCVariable.set(), slice_viewers, sounding_viewers, TEMPORARY_STRING_LENGTH, verbPrintf(), view_diffusivity, view_log_nuH, and view_nuH.

Referenced by misc_setup().

PetscErrorCode initBasalTillModel ( ) [protected, virtual]

Initialize the pseudo-plastic till mechanical model.

See IceBasalResistancePlasticLaw and updateYieldStressUsingBasalWater() and getBasalWaterPressure() and diffuseHmelt() for important model equations.

Calls either invertSurfaceVelocities(), for one way to get a map of till friction angle vtillphi, or computePhiFromBedElevation() for another way, or leaves vtillphi unchanged. First two of these are according to options -surf_vel_to_phi and -topg_to_phi, respectively.

Definition at line 53 of file iMbasal.cc.

References basal, IceGrid.com, computePhiFromBedElevation(), config, pism_config_editor.filename, NCConfigVariable.get(), NCConfigVariable.get_flag(), grid, PISMOptionsIsSet(), PISMOptionsString(), IceBasalResistancePlasticLaw.printInfo(), secpera, IceModelVec.set(), verbPrintf(), and vtauc.

Referenced by misc_setup(), and model_state_setup().

PetscErrorCode initFromFile ( const char * filename ) [virtual]

Read a saved PISM model state in NetCDF format, for complete initialization of an evolution or diagnostic run.

Before this is run, the method IceModel.grid_setup() determines the number of grid points (Mx,My,Mz,Mbz) and the dimensions (Lx,Ly,Lz) of the computational box from the same input file.

Reimplemented in IcePSTexModel, and IceMISMIPModel.

Definition at line 206 of file iMIO.cc.

References NCTool.close(), IceGrid.com, config, Enthnew3, NCTool.find_variable(), NCTool.get_att_text(), NCTool.get_dim_length(), NCConfigVariable.get_flag(), PISMVars.get_variables(), global_attributes, grid, have_ssa_velocities, vnreport.nc, NCTool.open_for_reading(), PISMOptionsIsSet(), NCGlobalAttributes.prepend_history(), IceGrid.rank, IceModelVec.read(), IceModelVec2V.read(), IceModelVec.set(), IceModelVec.set_attrs(), IceModelVec.set_name(), setEnth3FromT3_ColdIce(), setEnth3FromT3AndLiqfrac3(), tau3, variables, vel_ssa, and verbPrintf().

Referenced by model_state_setup().

PetscErrorCode initSSA ( ) [protected, virtual]

Each step of SSA uses previously saved values to start iteration; zero them here to start.

Definition at line 74 of file iMssa.cc.

References have_ssa_velocities, IceModelVec.set(), and vel_ssa.

Referenced by IceExactSSAModel.run(), and velocity().

bool issounding	(	const PetscInt	i,
		const PetscInt	j
	)			`[protected, virtual]`

Definition at line 357 of file iMutil.cc.

References jd.

Referenced by enthalpyAndDrainageStep().

PetscErrorCode massContExplicitStep ( ) [protected, virtual]

Update the thickness from the horizontal velocity and the surface and basal mass balance.

The partial differential equation describing the conservation of mass in the map-plane (parallel to the geoid) is

$\frac{\partial H}{\partial t} = M - S - \nabla\cdot \mathbf{q}$

where

$\mathbf{q} = \bar{\mathbf{U}} H.$

In these equations $H$ is the ice thickness, $M$ is the surface mass balance (accumulation or ablation), $S$ is the basal mass balance (e.g. basal melt or freeze-on), and $\bar{\mathbf{U}}$ is the vertically-averaged horizontal velocity of the ice. This procedure uses conservation of mass to update the ice thickness.

The PISMSurfaceModel pointed to by IceModel.surface provides $M$ . The PISMOceanModel pointed to by IceModel.ocean provides $S$ at locations below floating ice (ice shelves).

The map-plane flux of the ice $\mathbf{q}$ is defined by the above formula. Nonetheless the mass flux is split into the parts caused by non-sliding SIA-type deformation and caused by a nonzero basal sliding velocity:

$\mathbf{q} = - D \nabla h + \mathbf{U}_b H.$

Here $D$ is the (positive, scalar) effective diffusivity of the SIA and $\mathbf{U}_b$ is the basal sliding velocity.

The methods used are first-order explicit in time. The derivatives in $\nabla \cdot \mathbf{q}$ are computed by centered finite difference methods. In the case of the SIA contribution, the value of $D \nabla h$ is already stored in IceModelVec2Stag uvbar on the staggered grid by velocitySIAStaggered(). It is differenced in the standard centered manner (with averaging of the thickness onto the staggered grid).

Basal sliding may come from SSA or from a sliding law in SIA (the latter is usually inferior as a physical model). The divergence of $\mathbf{U}_b H$ is computed by upwinding after expanding

$\nabla\cdot (\mathbf{U}_b H) = \mathbf{U}_B \cdot \nabla H + (\nabla \cdot \mathbf{U}_B) H.$

That is, in the case of pure basal sliding the mass conservation equation is regarded as an advection equation with source term,

$\frac{\partial H}{\partial t} + \mathbf{U}_b \cdot \nabla H = M - S - (\nabla \cdot \mathbf{U}_b) H.$

The product of velocity and the gradient of thickness on the left is computed by first-order upwinding. Note that the CFL condition for this advection scheme is checked; see computeMax2DSlidingSpeed() and determineTimeStep().

Note that if the point is flagged as FLOATING_OCEAN0 then the thickness is set to zero. Note that the rate of thickness change $\partial H/\partial t$ is computed and saved, as is the rate of volume loss or gain.

Definition at line 328 of file iMgeometry.cc.

References acab, IceModelVec.add(), IceModelVec2Stag.begin_access(), IceModelVec.begin_access(), IceModelVec.beginGhostComm(), cell_area, check_maximum_thickness(), IceGrid.com, compute_ice_area(), computeSIAVelocities, config, IceModelVec.copy_to(), dt, dvoldt, IceGrid.dx, IceGrid.dy, IceModelVec2Stag.end_access(), IceModelVec.end_access(), IceModelVec.endGhostComm(), gdHdtav, IceModelVec2S.get_array(), NCConfigVariable.get_flag(), grid, PISMSurfaceModel.ice_surface_mass_flux(), IceModelVec2Mask.is_floating(), MASK_DRAGGING_SHEET, MASK_OCEAN_AT_TIME_0, ocean, pi, IceModelVec.scale(), secpera, PISMOceanModel.shelf_base_mass_flux(), shelfbmassflux, IceModelVec2S.sum(), surface, uvbar, IceModelVec.v, IceModelVec2Mask.value(), vbmr, vdHdt, vel_basal, vel_ssa, vH, vMask, vWork2d, IceModelVec.was_created(), IceGrid.xm, IceGrid.xs, IceGrid.year, IceGrid.ym, and IceGrid.ys.

Referenced by step().

PetscErrorCode misc_setup ( ) [virtual]

Miscellaneous initialization tasks plus tasks that need the fields that can come from regridding.

Reimplemented in IceMISMIPModel, and IceExactSSAModel.

Definition at line 538 of file iMinit.cc.

References IceGrid.com, config, correct_cell_areas(), NCConfigVariable.get_flag(), global_attributes, grid, init_extras(), init_snapshots(), init_timeseries(), init_viewers(), initBasalTillModel(), output_vars, NCConfigVariable.set_flag(), set_output_size(), and verbPrintf().

Referenced by init().

PetscErrorCode model_state_setup ( ) [virtual]

Sets the starting values of model state variables.

There are two cases:

1) Initializing from a PISM output file.

2) Setting the values using command-line options only (verification and simplified geometry runs, for example) or from a bootstrapping file, using heuristics to fill in missing and 3D fields.

Calls IceModel.regrid().

This function is called after all the memory allocation is done and all the physical parameters are set.

Calling this method should be all one needs to set model state variables. Please avoid modifying them in other parts of the initialization sequence.

Also, please avoid operations that would make it unsafe to call this more than once (memory allocation is one example).

Reimplemented in IceROSSModel, and IceUnitModel.

Definition at line 324 of file iMinit.cc.

References beddef, pism_config_editor.filename, grid, PISMBedDef.init(), initBasalTillModel(), initFromFile(), last_bed_def_update, PISMOptionsString(), regrid(), set_vars_from_options(), stampHistoryCommand(), updateSurfaceElevationAndMask(), variables, and IceGrid.year.

Referenced by init(), and run().

PetscErrorCode putTempAtDepth ( ) [protected, virtual]

Create a temperature field within ice and bedrock from given surface temperature and geothermal flux maps.

In bootstrapping we need to guess about the temperature within the ice and bedrock if surface temperature and geothermal flux maps are given. This rule is heuristic but seems to work well anyway. Full bootstrapping will start from the temperature computed by this procedure and then run for a long time (e.g. $10^5$ years), with fixed geometry, to get closer to thermomechanically coupled equilibrium. See the part of the User's Manual on EISMINT-Greenland.

Consider a horizontal grid point i,j. Suppose the surface temperature $T_s$ and the geothermal flux $g$ are given at that grid point. Within the corresponding column, denote the temperature by $T(z)$ for some elevation $z$ above the base of the ice. (Note ice corresponds to $z>0$ while bedrock has $z<0$ .) Apply the rule that $T(z)=T_s$ is $z$ is above the top of the ice (at $z=H$ ).

Within the ice, set

$T(z) = T_s + \alpha (H-z)^2 + \beta (H-z)^4$

where $\alpha,\beta$ are chosen so that

$\frac{\partial T}{\partial z}\Big|_{z=0} = - \frac{g}{k_i}$

and

$\frac{\partial T}{\partial z}\Big|_{z=H/4} = - \frac{g}{2 k_i}.$

The point of the second condition is our observation that, in observed ice, the rate of decrease in ice temperature with elevation is significantly decreased at only one quarter of the ice thickness above the base.

The temperature within the ice is not allowed to exceed the pressure-melting temperature.

Note that the above heuristic rule for ice determines $T(0)$ . Within the bedrock our rule is that the rate of change with depth is exactly the geothermal flux:

$T(z) = T(0) - \frac{g}{k_r} z.$

Note that $z$ here is negative, so the temperature increases as one goes down into the bed.

Definition at line 371 of file iMbootstrap.cc.

References artm, IceModelVec.begin_access(), IceModelVec.beginGhostComm(), IceFlowLaw.beta_CC_grad, bootstrapSetBedrockColumnTemp(), config, EC, IceModelVec.end_access(), IceModelVec.endGhostComm(), Enth3, g, NCConfigVariable.get(), IceModelVec2S.get_array(), NCConfigVariable.get_flag(), EnthalpyConverter.getEnthPermissive(), EnthalpyConverter.getPressureFromDepth(), grid, ice, PISMSurfaceModel.ice_surface_temperature(), IceFlowLaw.k, IceGrid.kBelowHeight(), IceFlowLaw.meltingTemp, IceGrid.Mz, pism_config_editor.result, IceFlowLaw.rho, setEnth3FromT3_ColdIce(), IceModelVec3.setInternalColumn(), surface, T3, Tb3, vbed, vGhf, vH, IceGrid.xm, IceGrid.xs, IceGrid.year, IceGrid.ym, IceGrid.ys, and IceGrid.zlevels.

Referenced by bootstrapFromFile(), and IceEISModel.set_vars_from_options().

PetscErrorCode readShelfStreamBCFromFile ( const char * filename ) [virtual]

Read certain boundary conditions from a NetCDF file, for diagnostic SSA calculations.

This is not really a bootstrap procedure, but it has to go somewhere.

For now it is only called using "pross".

Definition at line 155 of file iMbootstrap.cc.

References MaskInterp.allowed_levels, IceModelVec2Stag.begin_access(), IceModelVec.begin_access(), NCTool.close(), IceGrid.com, IceModelVec2V.create(), IceModelVec2S.create(), IceModelVec2Stag.end_access(), IceModelVec.end_access(), NCTool.find_variable(), g, IceModelVec2Mask.get_array(), IceModelVec2S.get_array(), PISMIO.get_grid_info_2d(), grid, IceModelVec.interpolation_mask, fill_missing.mask, MASK_FLOATING, MASK_SHEET, vnreport.nc, MaskInterp.number_allowed, NCTool.open_for_reading(), IceModelVec2V.regrid(), IceModelVec.regrid(), IceModelVec.set(), IceModelVec.set_attrs(), IceModelVec.set_glaciological_units(), update_viewers(), IceModelVec.use_interpolation_mask, uvbar, IceModelVec.v, verbPrintf(), vMask, IceModelVec.write_in_glaciological_units, IceGrid.xm, IceGrid.xs, IceGrid.ym, and IceGrid.ys.

Referenced by IceROSSModel.model_state_setup().

PetscErrorCode regrid ( ) [protected, virtual]

Manage regridding based on user options. Call IceModelVec.regrid() to do each selected variable.

For each variable selected by option -regrid_vars, we regrid it onto the current grid from the NetCDF file specified by -regrid_from.

The default, if -regrid_vars is not given, is to regrid the 3 dimensional quantities tau3, Tb3 and either T3 or Enth3. This is consistent with one standard purpose of regridding, which is to stick with current geometry through the downscaling procedure. Most of the time the user should carefully specify which variables to regrid.

Definition at line 312 of file iMIO.cc.

References NCTool.close(), IceGrid.com, config, pism_config_editor.filename, g, PISMVars.get(), NCConfigVariable.get(), NCConfigVariable.get_flag(), PISMIO.get_grid_info(), NCTool.get_vertical_dims(), grid, GRID_3D, GRID_3D_BEDROCK, IceModelVec.grid_type(), vnreport.nc, LocalInterpCtx.no_regrid_bedrock, LocalInterpCtx.no_regrid_ice, NCTool.open_for_reading(), PISMOptionsString(), IceModelVec.regrid(), IceModelVec.string_attr(), variables, verbPrintf(), grid_info.z_len, and grid_info.zb_len.

Referenced by model_state_setup().

PetscErrorCode report_grid_parameters ( ) [virtual]

Definition at line 287 of file iMutil.cc.

References IceGrid.bed_vertical_spacing, IceGrid.com, IceGrid.dx, IceGrid.dy, IceGrid.dz_fine, IceGrid.dzbMAX, IceGrid.dzbMIN, IceGrid.dzMAX, IceGrid.dzMIN, IceGrid.end_year, EQUAL, grid, IceGrid.ice_vertical_spacing, IceGrid.Lbz, IceGrid.Lx, IceGrid.Ly, IceGrid.Lz, IceGrid.Mbz, IceGrid.Mbz_fine, IceGrid.Mz, IceGrid.Mz_fine, IceGrid.printInfo(), IceGrid.printVertLevels(), IceGrid.start_year, verbPrintf(), and IceGrid.zblevels.

Referenced by init().

PetscErrorCode run ( ) [virtual]

Do the time-stepping for an evolution run.

This procedure is the main time-stepping loop.

Reimplemented in IceUnitModel, and IceExactSSAModel.

Definition at line 588 of file iceModel.cc.

References adaptReasonFlag, beddefEVENT, CFLmaxdt, CFLmaxdt2D, IceGrid.com, config, dt, dt_force, dt_from_cfl, dt_from_diffus, dtTempAge, dvoldt, IceGrid.end_year, endOfTimeStepHook(), gdHdtav, gDmax, NCConfigVariable.get_flag(), getVerbosityLevel(), global_attributes, gmaxu, gmaxv, gmaxw, grid, massbalEVENT, maxdt_temporary, model_state_setup(), PismLogEventRegister, PISMOptionsInt(), reportPATemps, secpera, NCVariable.set_string(), setVerbosityLevel(), siaEVENT, skipCountDown, ssaEVENT, IceGrid.start_year, stdout_flags, step(), summary(), summaryPrintLine(), tempEVENT, total_basal_ice_flux, total_sub_shelf_ice_flux, total_surface_ice_flux, update_viewers(), velmiscEVENT, verbPrintf(), write_extras(), write_snapshot(), write_timeseries(), and IceGrid.year.

Referenced by main().

PetscErrorCode set_grid_defaults ( ) [virtual]

Set default values of grid parameters.

Derived classes (IceCompModel, for example) reimplement this to change the grid initialization when no -i option is set.

Reimplemented in IceEISModel, IceMISMIPModel, IceUnitModel, IceCompModel, and IceExactSSAModel.

Definition at line 31 of file iMinit.cc.

References NCTool.close(), IceGrid.com, pism_config_editor.filename, NCTool.find_dimension(), NCTool.find_variable(), PISMIO.get_grid_info(), grid, grid_info.Lx, IceGrid.Lx, grid_info.Ly, IceGrid.Ly, mapping, vnreport.nc, NCTool.open_for_reading(), PISMOptionsIsSet(), PISMOptionsString(), NCConfigVariable.print(), NCConfigVariable.read(), secpera, grid_info.time, verbPrintf(), grid_info.x0, IceGrid.x0, grid_info.y0, IceGrid.y0, and IceGrid.year.

Referenced by grid_setup().

PetscErrorCode set_grid_from_options ( ) [virtual]

Initalizes the grid from options.

Reads all of -Mx, -My, -Mz, -Mbz, -Lx, -Ly, -Lz, -Lbz, -z_spacing and -zb_spacing. Sets corresponding grid parameters.

Reimplemented in IceMISMIPModel.

Definition at line 134 of file iMinit.cc.

References IceGrid.bed_vertical_spacing, IceGrid.com, IceGrid.compute_horizontal_spacing(), IceGrid.compute_vertical_levels(), EQUAL, grid, IceGrid.ice_vertical_spacing, IceGrid.Lbz, IceGrid.Lx, IceGrid.Ly, IceGrid.Lz, IceGrid.Mbz, IceGrid.Mx, IceGrid.My, IceGrid.Mz, PISMOptionsInt(), PISMOptionsList(), PISMOptionsReal(), and QUADRATIC.

Referenced by grid_setup().

PetscErrorCode set_output_size	(	string	option,
		string	description,
		string	default_value,
		set< string > &	result
	)			`[virtual]`

Assembles a list of variables corresponding to an output file size.

Definition at line 284 of file iMoptions.cc.

References IceGrid.com, config, NCConfigVariable.get_flag(), NCConfigVariable.get_string(), PISMVars.get_variables(), grid, PISMOptionsList(), and variables.

Referenced by init_snapshots(), and misc_setup().

PetscErrorCode set_time_from_options ( ) [protected, virtual]

Determine the run length, starting and ending years using command-line options.

Reimplemented in IceMISMIPModel.

Definition at line 570 of file iMinit.cc.

References IceGrid.com, config, IceGrid.end_year, NCConfigVariable.get(), grid, PISMOptionsReal(), IceGrid.start_year, and IceGrid.year.

Referenced by grid_setup().

PetscErrorCode set_vars_from_options ( ) [virtual]

Sets starting values of model state variables using command-line options.

Sets starting values of model state variables using command-line options and (possibly) a bootstrapping file.

In the base class there is only one case: bootstrapping.

Reimplemented in IceEISModel, IcePSTexModel, IceROSSModel, IceMISMIPModel, IceUnitModel, IceCompModel, and IceExactSSAModel.

Definition at line 371 of file iMinit.cc.

References bootstrapFromFile(), IceGrid.com, pism_config_editor.filename, grid, PISMOptionsString(), and verbPrintf().

Referenced by model_state_setup().

void setConstantNuHForSSA ( PetscScalar nuH ) [protected, virtual]

Definition at line 426 of file iceModel.cc.

References config, NCConfigVariable.set_flag(), SSAStrengthExtension.set_notional_strength(), and ssaStrengthExtend.

Referenced by setFromOptions().

PetscErrorCode setCTSFromEnthalpy ( IceModelVec3 & useForCTS ) [protected, virtual]

Compute the CTS field, CTS = E/E_s(p), from Enth3, and put in a global IceModelVec3 provided by user.

The actual cold-temperate transition surface (CTS) is the level set CTS = 1.

Does not communicate ghosts for IceModelVec3 useForCTS.

Definition at line 189 of file iMenthalpy.cc.

References IceModelVec.begin_access(), EC, IceModelVec.end_access(), Enth3, EnthalpyConverter.getCTS(), IceModelVec3.getInternalColumn(), EnthalpyConverter.getPressureFromDepth(), grid, IceGrid.Mz, IceModelVec.set_attrs(), IceModelVec.set_name(), vH, IceGrid.xm, IceGrid.xs, IceGrid.ym, IceGrid.ys, and IceGrid.zlevels.

Referenced by compute_cts().

PetscErrorCode setDefaults ( ) [protected]

Assigns default values to the many parameters and flags in IceModel.

The order of precedence for setting parameters in PISM is:

default values: Reasonable values to set up the model with are given in setDefaults() and in file pism/src/base/iMdefaults. setDefaults() is called in the constructor for IceModel. It would be reasonable to have setDefaults() read the defaults from a (default!) NetCDF file of a format so that others could be substituted.
derived class overrides: The constructor of a derived class can choose its own defaults for data members of IceModel (and its own data members). These will override the above.
command line options: The driver calls IceModel.setFromOptions() after the instance of IceModel (or a derived class) is constructed. setFromOptions() is virtual but should usually be called first if a derived class has a setFromOptions.

The input file (-i or -boot_from) will not contain (in Feb 2008 version of PISM) any values for the quantities which are set in setDefaults(). (There are parameters which can be set at the command line or by the input file, like grid.Mx. For -i the data file has the final word but for -boot_from the command line options have the final word.)

The defaults should be reasonable values under all circumstances or they should indicate missing values in some manner.

Definition at line 45 of file iMdefaults.cc.

References IceGrid.com, computeSIAVelocities, config, executable_short_name, NCConfigVariable.get(), global_attributes, gmaxu, gmaxv, gmaxw, grid, holdTillYieldStress, jd, leaveNuHAloneSSA, IceGrid.Mx, IceGrid.My, PISM_Revision, realAgeForGrainSize, reportPATemps, NCConfigVariable.set_flag(), NCVariable.set_string(), shelvesDragToo, ssaMatlabFilePrefix, ssaSystemToASCIIMatlab, standard_gravity, updateHmelt, useConstantTillPhi, and verbPrintf().

Referenced by IceModel().

PetscErrorCode setEnth3FromT3_ColdIce ( ) [protected, virtual]

Compute Enth3 from temperature T3 by assuming the ice has zero liquid fraction.

First this method makes sure the temperatures is at most the pressure-melting value, before computing the enthalpy for that temperature, using zero liquid fraction.

Because of how EnthalpyConverter.getPressureFromDepth() works, the energy content in the air is set to the value that ice would have if it a chunk of it occupied the air; the atmosphere actually has much lower energy content. It is done this way for regularity (i.e. dEnth/dz computations).

Because Enth3 gets set, does ghost communication to finish.

Definition at line 41 of file iMenthalpy.cc.

References IceModelVec.begin_access(), IceModelVec.beginGhostComm(), EC, IceModelVec.end_access(), IceModelVec.endGhostComm(), Enth3, EnthalpyConverter.getEnthPermissive(), IceModelVec3.getInternalColumn(), EnthalpyConverter.getPressureFromDepth(), grid, IceGrid.Mz, T3, vH, IceGrid.xm, IceGrid.xs, IceGrid.ym, IceGrid.ys, and IceGrid.zlevels.

Referenced by energyStep(), initFromFile(), putTempAtDepth(), IceExactSSAModel.set_vars_from_options(), IceCompModel.set_vars_from_options(), and IceMISMIPModel.set_vars_from_options().

PetscErrorCode setEnth3FromT3AndLiqfrac3 ( IceModelVec3 & Liqfrac3 ) [protected, virtual]

Compute Enth3 from temperature T3 and liquid fraction.

Because Enth3 gets set, does ghost communication to finish.

Definition at line 75 of file iMenthalpy.cc.

References IceModelVec.begin_access(), IceModelVec.beginGhostComm(), EC, IceModelVec.end_access(), IceModelVec.endGhostComm(), Enth3, EnthalpyConverter.getEnthPermissive(), IceModelVec3.getInternalColumn(), EnthalpyConverter.getPressureFromDepth(), grid, IceGrid.Mz, T3, vH, IceGrid.xm, IceGrid.xs, IceGrid.ym, IceGrid.ys, and IceGrid.zlevels.

Referenced by initFromFile().

PetscErrorCode setExecName ( const char * my_executable_short_name ) [virtual]

Definition at line 432 of file iceModel.cc.

References executable_short_name.

Referenced by main().

PetscErrorCode setFromOptions ( ) [virtual]

Read runtime (command line) options and alter the corresponding parameters or flags as appropriate.

A critical principle of this procedure is that it will not alter IceModel parameters and flags unless the user sets an option to do so. Thus this base class setFromOptions() can be called by a derived class after the derived class has set its own defaults for base class parameters and flags.

Also, this procedure should not allocate memory or create new objects using the new operator.

In fact this procedure only reads the majority of the options. Some are read in initFromOptions(), writeFiles(), and setStartRunEndYearsFromOptions(), among other places.

Note there are no options to directly set dx, dy, dz, Lbz, and year as the user should not directly set these grid parameters. There are, however, options for directly setting Mx, My, Mz, Mbz and also Lx, Ly, Lz.

Note that additional options are read by PISM{Atmosphere|Surface|Ocean}Model instances, including -pdd... and -d?forcing options.

Reimplemented in IceEISModel, IcePSTexModel, IceMISMIPModel, IceCompModel, and IceExactSSAModel.

Definition at line 46 of file iMoptions.cc.

References check_old_option_and_stop(), IceGrid.com, config, NCConfigVariable.flag_from_option(), NCConfigVariable.get_flag(), grid, holdTillYieldStress, ICE_HYBRID, iceFactory, jd, PISMOptionsInt(), PISMOptionsIsSet(), PISMOptionsReal(), PISMOptionsString(), realAgeForGrainSize, NCConfigVariable.scalar_from_option(), secpera, NCConfigVariable.set_flag(), setConstantNuHForSSA(), IceFlowLawFactory.setType(), shelvesDragToo, ssaMatlabFilePrefix, ssaSystemToASCIIMatlab, and verbPrintf().

Referenced by init().

PetscErrorCode setLiquidFracFromEnthalpy ( IceModelVec3 & useForLiquidFrac ) [protected, virtual]

Compute the liquid fraction corresponding to Enth3, and put in a global IceModelVec3 provided by user.

Does not communicate ghosts for IceModelVec3 useForLiquidFrac.

Definition at line 153 of file iMenthalpy.cc.

References IceModelVec.begin_access(), EC, IceModelVec.end_access(), Enth3, IceModelVec3.getInternalColumn(), EnthalpyConverter.getPressureFromDepth(), EnthalpyConverter.getWaterFraction(), grid, IceGrid.Mz, IceModelVec.set_attrs(), IceModelVec.set_name(), vH, IceGrid.xm, IceGrid.xs, IceGrid.ym, IceGrid.ys, and IceGrid.zlevels.

Referenced by compute_liqfrac().

PetscErrorCode setMaskSurfaceElevation_bootstrap ( ) [protected, virtual]

Determine surface and mask according to information in bootstrap file and options.

grid.year has to be valid at the time of this call.

Definition at line 268 of file iMbootstrap.cc.

References IceModelVec.beginGhostComm(), IceGrid.com, config, IceModelVec.end_access(), IceModelVec.endGhostComm(), NCConfigVariable.get(), IceModelVec2Mask.get_array(), IceModelVec2S.get_array(), NCConfigVariable.get_flag(), grid, ice, fill_missing.mask, MASK_FLOATING, MASK_OCEAN_AT_TIME_0, MASK_SHEET, ocean, IceFlowLaw.rho, PISMOceanModel.sea_level_elevation(), vbed, verbPrintf(), vH, vh, vMask, IceGrid.xm, IceGrid.xs, IceGrid.year, IceGrid.ym, and IceGrid.ys.

Referenced by bootstrapFromFile().

PetscErrorCode SigmaSIAToRegular ( ) [protected, virtual]

Put the volume strain heating (dissipation heating) onto the regular grid.

At the end of velocitySIAStaggered() the volume strain-heating $\Sigma$ is available on the staggered grid. This procedure averages it onto the regular grid. $\Sigma$ is used in the temperature equation.

Communication of ghosted values of Vec vSigma must occur between velocitySIAStaggered() and this procedure for the averaging to work.

Definition at line 543 of file iMsia.cc.

References IceModelVec.begin_access(), IceModelVec.end_access(), IceModelVec2S.get_array(), IceModelVec3.getInternalColumn(), grid, IceGrid.kBelowHeight(), IceGrid.Mz, IceModelVec3.setColumn(), Sigma3, Sigmastag3, vH, IceGrid.xm, IceGrid.xs, IceGrid.ym, and IceGrid.ys.

Referenced by velocity().

PetscErrorCode stampHistory ( string str ) [protected, virtual]

Get time and user/host name and add it to the given string.

Definition at line 144 of file iMutil.cc.

References global_attributes, pism_username_prefix(), and NCGlobalAttributes.prepend_history().

Referenced by endOfTimeStepHook(), stampHistoryCommand(), stampHistoryEnd(), and IceMISMIPModel.writeMISMIPFinalFiles().

PetscErrorCode stampHistoryCommand ( ) [protected, virtual]

Build a history string from the command which invoked PISM.

Definition at line 97 of file iMutil.cc.

References global_attributes, grid, PISM_Revision, NCGlobalAttributes.prepend_history(), IceGrid.size, stampHistory(), and TEMPORARY_STRING_LENGTH.

Referenced by model_state_setup().

PetscErrorCode stampHistoryEnd ( ) [protected, virtual]

Build the particular history string associated to the end of a PISM run.

Definition at line 123 of file iMutil.cc.

References IceGrid.com, e, grid, stampHistory(), and TEMPORARY_STRING_LENGTH.

Referenced by writeFiles().

PetscErrorCode step	(	bool	do_mass_continuity,
		bool	do_energy,
		bool	do_age,
		bool	do_skip,
		bool	use_ssa_when_grounded
	)			`[virtual]`

The contents of the time-step.

During the time-step we perform the following actions:

call additionalAtStartTimestep() to let derived classes do more

determine the maximum time-step boundary models can take

apply the time-step restriction from the -extra_{times,file,vars} mechanism:

update the yield stress for the plastic till model (if appropriate)

update the velocity field; in some cases the whole three-dimensional field is updated and in some cases just the vertically-averaged horizontal velocity is updated; see velocity()

determine the time step according to a variety of stability criteria; see determineTimeStep()

update the age of the ice (if appropriate)

update the enthalpy (or temperature) field according to the conservation of energy model based (especially) on the new velocity field; see energyStep()

compute fluxes through ice boundaries

update the thickness of the ice according to the mass conservation model; see massContExplicitStep()

call additionalAtEndTimestep() to let derived classes do more

Definition at line 441 of file iceModel.cc.

References adaptReasonFlag, additionalAtEndTimestep(), additionalAtStartTimestep(), ageStep(), bed_def_step(), beddef, beddefEVENT, determineTimeStep(), diffuseHmelt(), dt, dtTempAge, energyStep(), extras_max_timestep(), grid, ice_mass_bookkeeping(), massbalEVENT, massContExplicitStep(), PISMComponent.max_timestep(), maxdt_temporary, ocean, secpera, IceModelVec.set(), skipCountDown, stdout_flags, surface, tempEVENT, updateHmelt, updateSurfaceElevationAndMask(), updateYieldStressUsingBasalWater(), vdHdt, velocity(), and IceGrid.year.

Referenced by run().

PetscErrorCode summary	(	bool	tempAndAge,
		bool	useHomoTemp
	)			`[protected, virtual]`

Reimplemented in IceCompModel.

Definition at line 248 of file iMreport.cc.

References ageStats(), CFLviolcount, IceGrid.com, computeFlowUbarStats(), computeMaxDiffusivity(), config, dt, dvoldt, IceModelVec.end_access(), energyStats(), gdHdtav, gDmax, IceModelVec2S.get_array(), NCConfigVariable.get_flag(), getVerbosityLevel(), gmaxu, gmaxv, gmaxw, grid, IceGrid.Mx, IceGrid.My, IceGrid.Mz, secpera, stdout_flags, summaryPrintLine(), vH, volumeArea(), IceGrid.xm, IceGrid.xs, IceGrid.year, IceGrid.ym, and IceGrid.ys.

Referenced by IceExactSSAModel.run(), and run().

PetscErrorCode summaryPrintLine	(	PetscTruth	printPrototype,
		bool	tempAndAge,
		PetscScalar	year,
		PetscScalar	delta_t,
		PetscScalar	volume,
		PetscScalar	area,
		PetscScalar	meltfrac,
		PetscScalar	H0,
		PetscScalar	T0
	)			`[protected, virtual]`

Print a line to stdout which summarizes the state of the modeled ice sheet at the end of the time step.

Generally, two lines are printed to stdout, the first starting with a space and the second starting with the character 'S' in the left-most column (column 1).

The first line shows flags for which processes executed, and the length of the time step (and/or substeps under option -skip). See IceModel.run() for meaning of these flags.

If IceModel.printPrototype is TRUE then the first line does not appear and the second line has alternate appearance. Specifically, different column 1 characters are printed:

'P' line gives names of the quantities reported in the 'S' line, the "prototype", while
'U' line gives units of these quantities. This column 1 convention allows automatic tools to read PISM stdout and produce time-series. The 'P' and 'U' lines are intended to appear once at the beginning of the run, while an 'S' line appears at every time step. This base class version gives a report based on the information included in the EISMINT II intercomparison of ice sheet models[EISMINT00].

Note that the inputs volume and area to this method are in m^3 and m^2, respectively. Thus all inputs to this method are in MKS except for year.

The resulting numbers on an 'S' line have the following meaning in this base class version:

ivol is the total ice sheet volume
iarea is the total area occupied by positive thickness ice
thick0 is the ice thickness at the center of the computational domain
temp0 is the ice basal temperature at the center of the computational domain The last two can be interpreted as "sanity checks", because they give information about a location which may or may not be "typical".

For more description and examples, see the PISM User's Manual. Derived classes of IceModel may redefine this method and print alternate information. Use of DiagnosticTimeseries may be superior, however.

Reimplemented in IceMISMIPModel.

Definition at line 390 of file iMreport.cc.

References IceGrid.com, config, NCConfigVariable.get(), NCConfigVariable.get_flag(), getVerbosityLevel(), grid, secpera, stdout_flags, stdout_ssa, and verbPrintf().

Referenced by IceExactSSAModel.run(), run(), and summary().

PetscErrorCode surfaceGradientSIA ( ) [protected, virtual]

Compute the surface gradient in advance of the SIA velocity computation.

There are two methods for computing the surface gradient.

The default method is to directly differentiate the surface elevation $h$ by the Mahaffy method Mahaffy .

The alternative method, using option -eta which sets use_eta_transformation = true, is to transform the thickness to something more regular and differentiate that. We get back to the gradient of the surface by applying the chain rule. In particular, as shown in CDDSV for the flat bed and $n=3$ case, if we define

$\eta = H^{(2n+2)/n}$

then $\eta$ is more regular near the margin than $H$ . So we compute the surface gradient by

$\nabla h = \frac{n}{(2n+2)} \eta^{(-n-2)/(2n+2)} \nabla \eta + \nabla b,$

recalling that $h = H + b$ . This method is only applied when $\eta > 0$ at a given point; otherwise $\nabla h = \nabla b$ .

In both cases we are computing the gradient by finite differences onto a staggered grid. In the method with $\eta$ we apply centered differences using (roughly) the same method for $\eta$ and $b$ that applies directly to the surface elevation $h$ in the alternative method.

The resulting surface gradient on the staggered grid is in four Vecs, vWork2d[k] for k=0,1,2,3; recall there are two staggered grid points per regular grid point. These Vecs are used in velocitySIAStaggered(), basalSlidingHeatingSIA(), and horizontalVelocitySIARegular().

Definition at line 53 of file iMsia.cc.

References IceModelVec.beginGhostComm(), config, IceGrid.dx, IceGrid.dy, IceModelVec.end_access(), IceModelVec.endGhostComm(), IceFlowLaw.exponent(), IceModelVec2S.get_array(), NCConfigVariable.get_flag(), grid, ice, n, vbed, vh, vH, vWork2d, IceGrid.xm, IceGrid.xs, IceGrid.ym, and IceGrid.ys.

Referenced by computeMaxDiffusivity(), and velocity().

PetscErrorCode temperatureStep	(	PetscScalar *	vertSacrCount,
		PetscScalar *	bulgeCount
	)			`[protected, virtual]`

Takes a semi-implicit time-step for the temperature equation.

This method should be kept because it is worth having alternative physics, and so that older results can be reproduced. In particular, this method is documented by papers [BBL,BBssasliding]. The new browser page BOMBPROOF, PISM's numerical scheme for conservation of energy essentially documents the cold-ice-BOMBPROOF method here, but the newer enthalpy-based method is slightly different and (we hope) a superior implementation of the conservation of energy principle.

The conservation of energy equation written in terms of temperature is

$\rho c_p(T) \frac{dT}{dt} = k \frac{\partial^2 T}{\partial z^2} + \Sigma,$

where $T(t,x,y,z)$ is the temperature of the ice. This equation is the shallow approximation of the full 3D conservation of energy. Note $dT/dt$ stands for the material derivative, so advection is included. Here $\rho$ is the density of ice, $c_p$ is its specific heat, and $k$ is its conductivity. Also $\Sigma$ is the volume strain heating (with SI units of $J/(\text{s} \text{m}^3) = \text{W}\,\text{m}^{-3}$ ).

We handle horizontal advection explicitly by first-order upwinding. We handle vertical advection implicitly by centered differencing when possible, and retreat to implicit first-order upwinding when necessary. There is a CFL condition for the horizontal explicit upwinding [MortonMayers]. We report any CFL violations, but they are designed to not occur.

The vertical conduction term is always handled implicitly (i.e. by backward Euler).

We work from the bottom of the column upward in building the system to solve (in the semi-implicit time-stepping scheme). The excess energy above pressure melting is converted to melt-water, and that a fraction of this melt water is transported to the ice base according to the scheme in excessToFromBasalMeltLayer().

The method uses equally-spaced calculation but the methods getValColumn(), setValColumn() interpolate back-and-forth from this equally-spaced calculational grid to the (usually) non-equally spaced storage grid.

An instance of tempSystemCtx is used to solve the tridiagonal system set-up here.

In this procedure four scalar fields are modified: vHmelt, vbmr, Tb3, and Tnew3. But vbmr and Tb3 will never need to communicate ghosted values (horizontal stencil neighbors). The ghosted values for T3 are updated from the values in Tnew3 in the communication done by energyStep(). There is a diffusion model for vHmelt in diffuseHmelt() which does communication for vHmelt.

The (older) scheme cold-ice-BOMBPROOF, implemented here, is very reliable, but there is still an extreme and rare fjord situation which causes trouble. For example, it occurs in one column of ice in one fjord perhaps only once in a 200ka simulation of the whole sheet, in my (ELB) experience modeling the Greenland ice sheet. It causes the discretized advection bulge to give temperatures below that of the coldest ice anywhere, a continuum impossibility. So as a final protection there is a "bulge limiter" which sets the temperature to the surface temperature of the column minus the bulge maximum (15 K) if it is below that level. The number of times this occurs is reported as a "BPbulge" percentage.

Reimplemented in IceCompModel.

Definition at line 165 of file iMtemp.cc.

References allowAboveMelting, artm, tempSystemCtx.bed_thermal_c_p, tempSystemCtx.bed_thermal_k, tempSystemCtx.bed_thermal_rho, IceModelVec.begin_access(), IceFlowLaw.beta_CC_grad, IceFlowLaw.c_p, checkThinNeigh(), IceGrid.com, config, columnSystemCtx.ddratio(), tempSystemCtx.dtTemp, dtTempAge, IceGrid.dx, tempSystemCtx.dx, IceGrid.dy, tempSystemCtx.dy, IceGrid.dz_fine, tempSystemCtx.dzbEQ, tempSystemCtx.dzEQ, IceModelVec.end_access(), excessToFromBasalMeltLayer(), NCConfigVariable.get(), IceModelVec2S.get_array(), IceModelVec3.getValColumn(), IceModelVec3Bedrock.getValColumnPL(), grid, ice, tempSystemCtx.ice_c_p, tempSystemCtx.ice_k, tempSystemCtx.ice_rho, PISMSurfaceModel.ice_surface_temperature(), tempSystemCtx.initAllColumns(), IceModelVec2Mask.is_floating(), jd, IceFlowLaw.k, IceGrid.Mbz_fine, IceFlowLaw.meltingTemp, IceGrid.Mz_fine, columnSystemCtx.norm1(), ocean, PISMOptionsIsSet(), IceGrid.rank, IceFlowLaw.rho, secpera, tempSystemCtx.setBasalBoundaryValuesThisColumn(), columnSystemCtx.setIndicesAndClearThisColumn(), tempSystemCtx.setSchemeParamsThisColumn(), tempSystemCtx.setSurfaceBoundaryValuesThisColumn(), IceModelVec3.setValColumnPL(), IceModelVec3Bedrock.setValColumnPL(), PISMOceanModel.shelf_base_mass_flux(), PISMOceanModel.shelf_base_temperature(), shelfbmassflux, shelfbtemp, tempSystemCtx.Sigma, Sigma3, tempSystemCtx.solveThisColumn(), surface, tempSystemCtx.T, T3, tempSystemCtx.T3, tempSystemCtx.Tb, Tb3, Tnew3, tempSystemCtx.u, u3, tempSystemCtx.v, v3, IceModelVec2Mask.value(), vbmr, verbPrintf(), vGhf, vH, vHmelt, columnSystemCtx.viewColumnValues(), columnSystemCtx.viewSystem(), vMask, vRb, tempSystemCtx.w, w3, x, IceGrid.xm, IceGrid.xs, IceGrid.year, IceGrid.ym, IceGrid.ys, and IceGrid.zlevels_fine.

Referenced by energyStep().

PetscErrorCode testConvergenceOfNu	(	IceModelVec2S	vNuH[2],
		IceModelVec2S	vNuHOld[2],
		PetscReal *	norm,
		PetscReal *	normChange
	)			`[protected, virtual]`

Compares saved to current product of vertically-averaged viscosity times height, the quantity denoted $\bar\nu H$ . This comparison is used to determine if the outer iteration is converged.

Verification and PST experiments suggest that an $L^1$ criterion for convergence is best. For verification there seems to be little difference, presumably because the solutions are smooth and the norms are roughly equivalent on a subspace of smooth functions. For PST, the $L^1$ criterion gives faster runs with essentially the same results. Presumably that is because rapid (temporal and spatial) variation in $\bar\nu H$ occurs at margins, occupying very few horizontal grid cells. For the significant (e.g.~in terms of flux) parts of the flow, it is o.k. to ignore a bit of bad behavior at these few places, and $L^1$ ignores it more than $L^2$ (much less $L^\infty$ , which might not work at all).