#include <PISMSurface.hh>

Inheritance diagram for PSForceThickness:

Public Member Functions
	PSForceThickness (IceGrid &g, const NCConfigVariable &conf, PISMSurfaceModel *input)
virtual	~PSForceThickness ()
PetscErrorCode	init (PISMVars &vars)
virtual void	attach_atmosphere_model (PISMAtmosphereModel *input)
virtual PetscErrorCode	ice_surface_mass_flux (IceModelVec2S &result)
virtual PetscErrorCode	ice_surface_temperature (IceModelVec2S &result)
	Does not modify ice surface temperature.
virtual PetscErrorCode	max_timestep (PetscReal t_years, PetscReal &dt_years)
virtual void	add_vars_to_output (string keyword, set< string > &result)
	Adds ftt_modified_acab to "big" output files.
virtual PetscErrorCode	define_variables (set< string > vars, const NCTool &nc, nc_type nctype)
virtual PetscErrorCode	write_variables (set< string > vars, string filename)
Protected Attributes
string	input_file
PetscReal	alpha
IceModelVec2S *	ice_thickness
	current ice thickness produced by IceModel.
IceModelVec2S	target_thickness
IceModelVec2S	ftt_mask
IceModelVec2S	ftt_modified_acab

Detailed Description

A class implementing a modified surface mass balance which forces ice thickness to a given target by the end of the run.

Definition at line 281 of file PISMSurface.hh.

Constructor & Destructor Documentation

PSForceThickness	(	IceGrid &	g,
		const NCConfigVariable &	conf,
		PISMSurfaceModel *	input
	)		`[inline]`

Definition at line 283 of file PISMSurface.hh.

References alpha, PISMComponent::config, and ice_thickness.

virtual ~PSForceThickness ( ) [inline, virtual]

Definition at line 290 of file PISMSurface.hh.

Member Function Documentation

void add_vars_to_output	(	string	keyword,
		set< string > &	result
	)		`[virtual]`

Adds ftt_modified_acab to "big" output files.

Reimplemented from Modifier< PISMSurfaceModel >.

Definition at line 917 of file PISMSurface.cc.

References PISMComponent::add_vars_to_output(), and Modifier< PISMSurfaceModel >::input_model.

void attach_atmosphere_model ( PISMAtmosphereModel * input ) [virtual]

Reimplemented from PSModifier.

Definition at line 656 of file PISMSurface.cc.

References PISMSurfaceModel::attach_atmosphere_model(), and Modifier< PISMSurfaceModel >::input_model.

PetscErrorCode define_variables	(	set< string >	,
		const NCTool &	,
		nc_type
	)		`[virtual]`

Defines requested couplings fields to file and/or asks an attached model to do so.

Reimplemented from Modifier< PISMSurfaceModel >.

Definition at line 926 of file PISMSurface.cc.

References PISMSurfaceModel::define_variables(), ftt_mask, ftt_modified_acab, Modifier< PISMSurfaceModel >::input_model, and target_thickness.

PetscErrorCode ice_surface_mass_flux ( IceModelVec2S & result ) [virtual]

If -force_to_thk foo.nc is in use then vthktarget will have a target ice thickness map. Let $H_{\text{tar}}$ be this target thickness, and let $H$ be the current model thickness. Recall that the mass continuity equation solved by IceModel::massContExplicitStep() is

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

and that this procedure is supposed to produce $M$ . In this context, the semantics of -force_to_thk are that $M$ is modified by a multiple of the difference between the target thickness and the current thickness. In particular, the $\Delta M$ that is produced here is

$\Delta M = \alpha (H_{\text{tar}} - H)$

where $\alpha>0$ is determined below. Note $\Delta M$ is positive in areas where $H_{\text{tar}} > H$ , so we are adding mass there, and we are ablating in the other case.

Let $t_s$ be the start time and $t_e$ the end time for the run. Without flow or basal mass balance, or any surface mass balance other than the $\Delta M$ computed here, we are solving

$\frac{\partial H}{\partial t} = \alpha (H_{\text{tar}} - H)$

Let's assume $H(t_s)=H_0$ . This initial value problem has solution $H(t) = H_{\text{tar}} + (H_0 - H_{\text{tar}}) e^{-\alpha (t-t_s)}$ and so

$H(t_e) = H_{\text{tar}} + (H_0 - H_{\text{tar}}) e^{-\alpha (t_e-t_s)}$

The constant $\alpha$ has a default value pism_config:force_to_thickness_alpha.

The final feature is that we turn on this mechanism so it is harshest near the end of the run. In particular,

$\Delta M = \lambda(t) \alpha (H_{\text{tar}} - H)$

where

$\lambda(t) = \frac{t-t_s}{t_e-t_s}$

The next exacmple uses files generated from the EISMINT-Greenland experiment; see the corresponding chapter of the User's Manual.

Suppose we regard the SSL2 run as a spin-up to reach a better temperature field. It is a spinup in which the surface was allowed to evolve. Assume the early file green20km_y1.nc has the target thickness, because it essentially has the input thickness. This script adds a 500 a run, to finalize the spinup, in which the ice sheet geometry goes from the the thickness values in green_ssl2_110ka.nc to values very close to those in green20km_y1.nc:

#!/bin/bash

NN=8

mpiexec -n $NN pismr -ys -500.0 -ye 0 -skip 5 -i green_ssl2_110ka.nc -atmosphere searise_greenland \
  -surface pdd,forcing -pdd_fausto -force_to_thk green20km_y1.nc \
  -o with_force.nc -ts_file ts_with_force.nc -ts_times -500:10:0
  
mpiexec -n $NN pismr -ys -500.0 -ye 0 -skip 5 -i green_ssl2_110ka.nc -atmosphere searise_greenland \
  -surface pdd -pdd_fausto \
  -o no_force.nc -ts_file ts_no_force.nc -ts_times -500:10:0

mpiexec -n $NN pismr -ys -500.0 -ye 0 -skip 5 -i green_ssl2_110ka.nc -atmosphere searise_greenland \
  -surface pdd,forcing -pdd_fausto -force_to_thk green20km_y1.nc -force_to_thk_alpha 0.0002 \
  -o weak_force.nc -ts_file ts_weak_force.nc -ts_times -500:10:0

The script also has a run with no forcing, and one with forcing at a lower alpha value, a factor of ten smaller than the default.

As shown below, the time series for ivol in the above time series files show that the force-to-thickness mechanism is forcing a system with negative feedback. We see decaying oscillations toward the intended volume.

Volume results from the -force_to_thk mechanism.

Reimplemented from PSModifier.

Definition at line 838 of file PISMSurface.cc.

References alpha, ftt_mask, ftt_modified_acab, PISMComponent::grid, exactV::H(), PISMSurfaceModel::ice_surface_mass_flux(), ice_thickness, Modifier< PISMSurfaceModel >::input_model, PISMComponent_TS::t, and target_thickness.

PetscErrorCode ice_surface_temperature ( IceModelVec2S & result ) [virtual]

Does not modify ice surface temperature.

Reimplemented from PSModifier.

Definition at line 882 of file PISMSurface.cc.

References PISMSurfaceModel::ice_surface_temperature(), and Modifier< PISMSurfaceModel >::input_model.

PetscErrorCode init ( PISMVars & vars ) [virtual]

Reimplemented from PISMSurfaceModel.

Definition at line 660 of file PISMSurface.cc.

References alpha, ftt_mask, ftt_modified_acab, PISMComponent::grid, ice_thickness, PISMSurfaceModel::init(), input_file, Modifier< PISMSurfaceModel >::input_model, vnreport::nc, secpera, and target_thickness.

PetscErrorCode max_timestep	(	PetscReal	t_years,
		PetscReal &	dt_years
	)		`[virtual]`

The timestep restriction is, by direct analogy, the same as for

$\frac{dy}{dt} = - \alpha y$

with explicit (forward) Euler. If $\Delta t$ is the time step then Euler is $y_{n+1} = (1-\alpha \Delta t) y_n$ . We require for stability that $|y_{n+1}|\le |y_n|$ , which is to say $|1-\alpha \Delta t|\le 1$ . Equivalently (since $\alpha \Delta t>0$ ),