Model Structure

skadia.f90

The main program calls readcontrol to read in the user’s preferences for the model run from the control file control.nml. The information in this file determines several factors within the model; the length of the model run, the mode the model is going to be run in, and the input and output files required and generated, respectively.

Call Readoptimfile reads in the initial and the range of values that can be taken for the parameters that are able to be optimized within the model. Depending on the mode that the model is running in determines which set of routines are called next,

1. Mode = 0; Optimisation mode
2. Mode = 1; Phase space mode
3. Mode = 2; Inverse model mode
4. Mode = 3; Forward mode

Valley

This is the main program that runs the valley glacier model. Within this routine the following subroutines and functions are called before the model time step is initiated:

1. call getoptimvalues
2. call initial
3. call readgeometry
4. function calcwdth
5. call readcostdata

Main time loop is controlled by count and determined variables trun, tinc and tstr (all units of years and defined in the control file).

Call getoptimvalues

Depending on which mode the model is running in the parameter values for the current model run are selected. In the optimisation mode, the final model run uses the 'best' parameter set. In the other modes, this routine selects parameter values from either; the parameter file, the space survey or the genetic algorithm (GaFortran).

Call initial

This routine is where the parameter values are scaled. For further information see the section, scaling.

Call readgeometry

to read basic geometry of glacier system. The format and name of the input file containing the geometry information is shown in the control file. The information to be read in:

• current surface and bedrock elevations
• reference widths
• reference thicknesses.

All the variables are scaled and the current ice thickness is determined from the upper surface and bedrock elevation.

Call readcostdata

The files that contain glacier observations against which the model's performance is measured are read in. The names and the format of these files are described in the control file section. The model is capable of using a range of optimisation data including;

• ice thickness profiles through time
• seasonal and annual/net mass balance for the whole glacier
• seasonal and annual/net mass balance for different elevation bands.
• glacier length changes
• ela changes.

---

<font-size: 20>Within the main time loop the following routines are called from the main Valley program.

Call readmet

Reads in the meteorological data for the current year. The format of the meteorological data read in depends on the mass balance method being used in the model. In its current configuration the model requires the following fields;

• year
• julian day
• temperature
• temperature range
• precipitation.
  

Each year's data is in a seperate text file.

Function calcacab

This function returns the mass balance profile for the current year. Depending on the setting of whichacab, the mass balance is either calculated using a degree day method (whichacab set to 0) or a simple energy balance model (whichacab set to 1). To account for superimposed ice a simple model is used whereby a fraction of the meltwater refreezes in the snowpack, after which susequent melt runs off. The routines that are used in this function are described in the section Model mass balance.

The function is called once per year from the main valley model and has a daily time step.

1. There are stores for the depths of snow, super-imposed ice and glacier ice along the glacier. The latter is set the zero every year and is a relative store, while the other two inherit depths from the end of the previous year.
2. The main loops are for days in the year (leap years are determined elsewhere and dayinyr is corrected for them) and points along the glacier. Mass balance is determined at end point on the transect (glacier or not).
3. The loop proceeds by determining surface air temperatures using a lapse rate and the daily mean temperature at the met. station (or equivalent model-based data or assumed elevation).
4. A similar procedure is used to estimate daily precipitation totals along the glacier. Note that this is only done if precipitation actually occurred on the day in question.
5. An estimation of the amount of precipitation falling as snow is then made based on diurnal air temperatures.
6. Total potential melt rates are then determined for the day using either day-degree or a simple energy-balance model. In the former case, snow is initially assumed but a correction for ice is made within the surface mass balance subroutine.
7. The evolution of the snow, superimposed ice and glacier ice along the transect is then determined using the potential melt totals, as well as daily rainfall and total precipitation.
8. After the main loops, we tidy up by remembering the end of year depths of snow and super-imposed ice and sorting out the mass-balance season information required in the cost function.

Function calcmask

This function is called from valley and returns a mask that indicates where calculations for ice flow etc need to be performed and where not. In the former case, a unique integer value is given to each cell; while 0 is entered in the latter case. For a point to be included one of the following must be true:

• the point has ice thickness greater than zero;
• its mass balance is positive;
• one of its up- or down-stream neighbours has non-zero ice thickness; or
• the point is at a branch in the system and one of its branched neighbours has non-zero ice thickness.

<font-size: 20>If thickness/width evolution is included (whichthck =1) then the model runs a set of rotuines to determine the glacier's dynamic response (Model dynamics).