In this exercise, you will use CESM to run a Greenland ice sheet simulation.

You will
modify the configuration file of the ice sheet model to add a surface mass balance
anomaly to the ice sheet which will induce ice loss particularly along the coast of the
Greenland ice sheet. This experiment is based on the initMIP Greenland experiment
which was a precursor to the ISMIP6 Greenland experiment (part of the CMIP6
contribution).

NOTE: you will be using an experimental version of CESM that has not been tested to run
experiments outside of the one proposed in this document.

Goal of the experiment

You will create a case that consists of an active ice-sheet component (CISM, the Community
Ice Sheet Model), and stub atmosphere, land, river transport, ocean, sea ice and wave
components. This case will show you an example on how to run a single component within
the CESM framework. The experiment begins with a spun-up Greenland Ice Sheet (GrIS) in
equilibrium with a Surface Mass Balance (SMB) climatology from RACMO 2.3 spanning the
years 1958-2015. You will then apply a SMB anomaly that will promote mass loss over
GrIS. This experiment is similar to the one performed for the initMIP Greenland experiment
which was a precursor to the ISMIP6 Greenland stand-alone experiment (which is one of
the ice sheet component simulation for the IPCC-AR6 report).

CESM version for this experiment

Typically, CISM requires at least to run with a data land model. Here, the surface mass
balance anomaly is already included in the input file needed to run the experiment. This
experiment could also be run using CISM stand alone, (i.e, outside of the CESM framework).
This version of CESM is called “cesm_cism_tutorial” and is not available on the CESM github
repository but is available in the directory “/glade/u/home/gunterl/CESM_tutorial_2019/”

Note: because this CESM version is experimental, some glitches still exist but will not
impact the model simulation. For instance, the surface mass balance and air temperature
fields are not being written properly to the output file. This has very little influence on the
simulation and its outcome.

Create a case

> cd /glade/u/home/gunterl/CESM_tutorial_2019/cesm_cism_tutorial/cime/scripts

> ./create_newcase --case ~/cases/initMIP_GrIS.1 --mach cheyenne --res

f19_g17 --compset 1850_SATM_SLND_SICE_SOCN_SROF_CISM2%EVOLVE_SWAV --run-
unsupported --project UESM0011
The resolution f19_g17 implicitly sets the CISM grid to its default of a 4-km Greenland grid
(if you’d like, you can use the equivalent alias, f19_g17_gl4, in order to be more explicit).

You need the ‘--run-unsupported’ option here because CESM doesn’t have any automated
tests of this exact compset and resolution combination. If you’d like, try the above
command without ‘--run-unsupported’ and see what happens. In general, you should be
wary of running configurations where you need to specify ‘--run-unsupported’, because
you’re more likely to run into problems (obvious or subtle), but in this case I’ll vouch for
this being an okay configuration.
Note: remember that the project code UESM0011 is only available for the tutorial. You will
need to replace this code by one that contains your HPC allocation in the future.

Some relevant information printed by create_newcase is:

- GLC component is cism2 (default, higher-order, can run in
parallel):cism ice evolution turned on:
o CISM2 gives us the latest suite of CISM dynamics and parameterizations.
o Grid: g%gland4 which indicates that the run will use the 4 km resolution GrIS
grid (which is what we want).

Configure and build the model

> cd ~/cases/initMIP_GrIS.1

Now setup and build the case.

> ./case.setup

Now we need to edit the “user_nl_cism” file to run with a specific input file and apply the
anomaly field to the run. The reason for this new input file is because we typically apply
external forcing using data model within CESM. The anomaly field that we will apply in this
experiment has been added to the default input file and we could access it easily via a
simple change in the CISM configuration file.

Open the “user_nl_cism” file:

> vim user_nl_cism

Add the following lines:

cisminputfile =
‘/glade/u/home/gunterl/CESM_tutorial_2019/input_initMIP_GrIS/greenland_4km_ep
sg3413_c171126_initMIP_anomaly.nc’
enable_acab_anomaly = .true.

WARNING: if you simply copy and paste the commands above, please make sure that you
retype the quotes after “cisminputfile” directly from your keypad in your terminal. MS
Word changes their shapes and cesm does not like that.
To make sure there were no typos, run

> ./preview_namelists

and make sure that command completes successfully.

Now run:
> qcmd -A UESM0011 -- ./case.build

Note: typically, you could replace the command above by simply typing “./case.build”.

When the code is done with the build phase with no errors you should see a message:
“MODEL BUILD HAS FINISHED SUCCESSFULLY”

Submit the run

It's often a good idea to first do a short test run, to make sure you haven't introduced any
errors that prevent the model from running. For this case let’s trust what we have done and
go ahead start a 100-year run. Note that we are only doing this here because we are
running a cheap CISM only run. For more involved experiments it is highly recommended
to perform a test run first.

Let’s setup a 100-year run:

> ./xmlchange STOP_OPTION=nyears,STOP_N=100

Also, lengthen the wall-clock limit to allow enough time for this run:

> ./xmlchange JOB_WALLCLOCK_TIME=0:40

Then submit the run:

> ./case.submit

You can check on the status of you run with:

> qstat -u $USER

Wait for the run to finish

The run will take 10-20 min. to complete. While you’re waiting, you can take some time to
look through these two things:

(1) The various namelist settings that can be used to control CISM’s evolution:
http://www.cesm.ucar.edu/models/cesm2/settings/current/cism_nml.html

The column “Variable” gives you the name of a configuration that is available in CISM and
can be written in the user_nl_cism similarly to what we just did. The “Namelist Group”
indicates the group in the namelist file under which the variable is written.
Note that the options under headings beginning with “cism_config” will appear in the
cism.config file in your run directory; options in the other sections (which do not mention
“cism_config”) will appear in cism_in. However, both sets of options are modified in the
same way: via the user_nl_cism file in your case directory.

(2) The CISM source code, as described in the Appendix.

View the results

Now we will look at some output from your run. When a run completes successfully,
output is written to several archive directories. Let's look at some glc output:

> cd /glade/scratch/$USER/archive/initMIP_GrIS.1/glc/hist
> ls

You should see one file for each year of the run. Each of the yearly files contains a number
of yearly average fields. This command will tell you about the contents of the first yearly
file:

> module load netcdf

> ncdump -h initMIP_GrIS.1.cism.initial_hist.0001-01-01-00000.nc| less

The file contains several output variables. Remember, in this experiment we added a
surface mass balance anomaly to be applied to a constant surface mass balance. In CISM,
the surface mass balance can be read as “smb” (surface mass balance) or “acab”
(accumulation, ablation rate). Both fields differ by their units. “acab” has units of m/yr ice
equivalent and “smb” has units of mm/yr water equivalent. This distinction comes from
that an ice sheet model typically computes ice related variables while an atmospheric
model computes precipitation in water content. The conversion from one to another can be
done by using ice and water densities.
While “smb” is an output variable, it returns 0 because of a glitch in this CESM version. The
output file also returns “acab_applied” which is similar to “acab” but is the actual applied
accumulation, ablation to the ice sheet (including the anomaly, a filter setting 0 surface
mass balance in the open ocean, …). Values for acab (or smb) are positive for ice
accumulation (precipitation) and negative for ice ablation (melting).

You can look at individual files using the ncview command. For example, to look at results
for year 50:

> module load ncview

> ncview initMIP_GrIS.1.cism.h.0050-01-01-00000.nc

Click on the box labeled “acab_applied” to look at the applied accumulation. You can also
click on “thk” to look at the actual thickness of the ice sheet at that year.

CISM outputs yearly scalar values in the history files:

imass: the mass of ice available (kg)
imass_above_flotation: the mass above flotation available (kg)
iareaf: area covered by floating mask (m2)
iareag: area covered by grounded mask (m2)
total_bmb_flux: total basal mass balance flux (kg/s)
total_calving_flux: total calving mass balance flux (kg/s)
total_gl_flux: total grounding line flux (kg/s)
total_smb_flux: total surface mass balance flux (kg/s)

As such it is not convenient to look at them. Let us combine them and look at their time
series. To do so we will use the NetCDF operator (also called NCO):

> module load nco

> ncrcat -v
imass,imass_above_flotation,iareaf,iareag,total_bmb_flux,total_calving_flux,t
otal_gl_flux,total_smb_flux initMIP_GrIS.1.cism.initial_hist.0001-01-01-
00000.nc initMIP_GrIS.1.cism.h.*.nc out_scalars.initMIP_GrIS.nc

This command will concatenate all of the scalar output from the files
“initMIP_GrIS.cism.initial_hist.0001-01-01-00000.nc” and “initMIP_GrIS.cism.h.*.nc” (where
* is the wildcard character which will look for all the files in the directory that have the
same pattern. In this specific case, it really replaces the numbers 0002 to 0101) and will
place them all in a file called “out_scalars.initMIP_GrIS.nc”. Now take a look at that file:

> ncview out_scalars.initMIP_GrIS.nc

If you click on “imass” you will see a time series that decreases over time indicating the ice
sheet is losing mass. If you look at “imass_above_flotation” you will see a similar curve as
for “imass”. This particular scalar is useful in diagnosing the contribution of the ice sheet to
sea level change – floating mass loss does not contribute to sea level change, only ice
change that is grounded and not already in the ocean (i.e above sea level). As an exercise
you could calculate the amount of sea level contribution from this experiment using the fact
that 1mm of eustatic sea level rise is equivalent to 360 Gt of ice (or 360.e12 kg). Your
answer should yield 160.63 mm contribution to sea level rise.

You could also concatenate all variables in the history files without specifying the variable
names:

> ncrcat initMIP_GrIS.1.cism.initial_hist.0001-01-01-00000.nc

initMIP_GrIS.1.cism.h.*.nc out_history.initMIP_GrIS.nc
(if you type this command you can ignore the warning about scale_factor.)

However, this will result in a very big file (about 11Gb). It might be more convenient to
move a single file around but it might not be convenient for visualization purposes as
ncview tends to be slower with larger files. Instead, I suggest you to look at the changes in
given fields between the end and the initial state of the ice sheet. You can create such a file
with:

> ncdiff initMIP_GrIS.1.cism.h.0101-01-01-00000.nc

initMIP_GrIS.1.cism.initial_hist.0001-01-01-00000.nc out.initMIP_GrIS.diff.nc

This command will create a file “out.initMIP_GrIS.diff.nc” which contains the difference
between the files “initMIP_GrIS.cism.h.0101-01-01-00000.nc” and
“initMIP_GrIS.cism.initial_hist.0001-01-01-00000.nc”. Take a look at it:

> ncview out.initMIP_GrIS.diff.nc

And click on the “thk” field. You should observe a decrease in ice thickness (negative
values) in most of the ice sheet domain especially around the ice sheet margins. You will
also notice a thickness increase in some areas in Eastern Greenland. Now click on the
“acab_applied” field. This difference actually displays the anomaly that has been applied
throughout the run with mostly negative values everywhere indicating this experiment was
meant to reduce ice sheet mass. If that is the case, is it surprising to see a positive ice
thickness difference anywhere on the ice sheet?

In order to answer this question fully we would first need to run a control experiment and
make sure that the ice sheet does not have drift in its system. Let’s assume it does not, then
I would invite you to select a region of interest and look at the ice thickness history. Ice
advection is a possible explanation. Upstream thicker ice moves downstream and replaces
thinner ice. Could you verify this assumption? (You might want to check ice velocities and
bed topography to get a full picture of what is happening).

Appendix: Ice sheet model directory structure

Go back to the CESM code:

> cd /glade/p/cesm/tutorial/cesm2_0_alpha07c/components
> ls

The components directory contains subdirectories corresponding to each possible active

CESM component.

Let's take a quick tour of “cism”:

> cd cism
> ls

The source code resides in several subdirectories:

- source_cism = The core CISM code; especially see the following subdirectories:
o libglissade: the glissade dynamical core, which is the default when running
with CISM2
o libglide: the serial, shallow ice approximation-based glide dynamical core;
this also contains some modules shared between glide and glissade
o libglimmer: lower-level, shared utilities
- source_glc = code needed to link CISM to the CESM coupler
- drivers = more code associated with the coupler

Feel free to explore these directories while waiting for the model to build or run.