◆ adjust_mass_flux()

void pism::surface::ForceThickness::adjust_mass_flux	(	double	time,
		const array::Scalar &	ice_thickness,
		const array::CellType &	cell_type,
		array::Scalar &	result
	)		const

private

If -force_to_thickness_file 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 GeometryEvolution 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_thickness 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:surface.force_to_thickness.alpha.

The next example 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  # default number of processors
if [ $# -gt 0 ] ; then  # if user says "test_ftt.sh 8" then NN = 8
  NN="$1"
fi
 
# set MPIDO if using different MPI execution command, for example:
#  $ export PISM_MPIDO="aprun -n "
if [ -n "${PISM_MPIDO:+1}" ] ; then  # check if env var is already set
  echo "$SCRIPTNAME      PISM_MPIDO = $PISM_MPIDO  (already set)"
else
  PISM_MPIDO="mpiexec -n "
  echo "$SCRIPTNAME      PISM_MPIDO = $PISM_MPIDO"
fi
 
# check if env var PISM_DO was set (i.e. PISM_DO=echo for a 'dry' run)
if [ -n "${PISM_DO:+1}" ] ; then  # check if env var DO is already set
  echo "$SCRIPTNAME         PISM_DO = $PISM_DO  (already set)"
else
  PISM_DO=""
fi
 
# prefix to pism (not to executables)
if [ -n "${PISM_PREFIX:+1}" ] ; then  # check if env var is already set
  echo "$SCRIPTNAME     PISM_PREFIX = $PISM_PREFIX  (already set)"
else
  PISM_PREFIX=""    # just a guess
  echo "$SCRIPTNAME     PISM_PREFIX = $PISM_PREFIX"
fi
 
# set PISM_EXEC if using different executables, for example:
#  $ export PISM_EXEC="pismr -energy cold"
if [ -n "${PISM_EXEC:+1}" ] ; then  # check if env var is already set
  echo "$SCRIPTNAME       PISM_EXEC = $PISM_EXEC  (already set)"
else
  PISM_EXEC="pismr"
  echo "$SCRIPTNAME       PISM_EXEC = $PISM_EXEC"
fi
 
 
PISM="${PISM_PREFIX}${PISM_EXEC}"
 
cmd="$PISM_MPIDO $NN $PISM -ys -1000.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 -1000:yearly:0"
$PISM_DO $cmd
 
echo
 
cmd="$PISM_MPIDO $NN $PISM -ys -1000.0 -ye 0 -skip 5 -i green_ssl2_110ka.nc -atmosphere searise_greenland \
  -surface pdd,forcing -pdd_fausto -force_to_thickness_file green20km_y1.nc \
  -o default_force.nc -ts_file ts_default_force.nc -ts_times -1000:yearly:0"
$PISM_DO $cmd
 
echo
 
cmd="$PISM_MPIDO $NN $PISM -ys -1000.0 -ye 0 -skip 5 -i green_ssl2_110ka.nc -atmosphere searise_greenland \
    -surface pdd,forcing -pdd_fausto -force_to_thickness_file green20km_y1.nc -force_to_thickness_alpha 0.005 \
    -o weak_force.nc -ts_file ts_weak_force.nc -ts_times -1000:yearly:0"
$PISM_DO $cmd
 
 
cmd="$PISM_MPIDO $NN $PISM -ys -1000.0 -ye 0 -skip 5 -i green_ssl2_110ka.nc -atmosphere searise_greenland \
    -surface pdd,forcing -pdd_fausto -force_to_thickness_file green20km_y1.nc -force_to_thickness_alpha 0.05 \
    -o strong_force.nc -ts_file ts_strong_force.nc -ts_times -1000:yearly:0"
$PISM_DO $cmd

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

Definition at line 219 of file ForceThickness.cc.

References pism::array::CellType::grounded(), m_alpha, m_alpha_ice_free_factor, pism::Component::m_config, m_ftt_mask, pism::Component::m_grid, m_ice_free_thickness_threshold, pism::Component::m_log, m_start_time, m_target_thickness, and pism::Component::time().

Referenced by update_impl().