Difference between revisions of "Model Structure"

Latest revision as of 17:12, 9 November 2007

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,

Mode = 0; Optimisation mode
Mode = 1; Phase space mode
Mode = 2; Inverse model mode
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:

call getoptimvalues
call initial
call readgeometry
function calcwdth
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.

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.
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).
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).
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.
An estimation of the amount of precipitation falling as snow is then made based on diurnal air temperatures.
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.
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.
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).

@@ Line 20: / Line 20: @@
-== Call getoptimvalues ==
+=== 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 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:
+=== 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
@@ Line 30: / Line 35: @@
 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.
-== 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.
-== function calcthck ==
-This function returns ice thickness along the glacier at the new time including the effects of local mass balance and ice flow.  This is done using an implicit technique and is complicated by the need to incorporate the effect of possible branches in the glacier system.  The equation to be solved is,
-::<math>\frac{\delta\,s} {\delta\,t}= b - \frac{1} {w}  \frac{\delta\,wq} {\delta\,x}</math> <br>
-:::<math>= b- \frac{\delta\,q} {\delta\,x}- \frac{q}{w} \frac{\delta\,w} {\delta\,x}</math> <br>
+<font-size: 20>'''Within the main time loop the following routines are called from the main Valley program.'''</font>
-:::<math>= b+ \frac{\delta} {\delta\,x}D \frac{\delta\,s} {\delta\,x}+ \frac{D}{w} \frac{\delta\,w} {\delta\,x}\frac{\delta\,s} {\delta\,x}</math> <br>
-:::<math>= b+ \frac{\delta} {\delta\,x}D \frac{\delta\,s} {\delta\,x}+ c \frac{\delta\,s} {\delta\,x}</math>. <br>
+=== 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.<br>
-Note that we are not incorporating changes in bedrock elevation so that <math>\frac{\delta\,s} {\delta\,t}= \frac{\delta\,H} {\delta\,t} </math>. The numerical technique involved is based on Crank-Nicolson time differencing:
+Each year's data is in a seperate text file.
-::<math> s_{t+\Delta\,t}=s_{t}+ \Delta\,t \left [{b}\right ]_{t}+ \Delta\,t\alpha\left [{\frac{\delta}{\delta\,x}D\frac{\delta\,s} {\delta\,x}+c\frac{\delta\,s} {\delta\,x}}\right ]_{t}+\Delta\,t \left({1-\alpha}\right)\left [{\frac{\delta}{\delta\,x}D\frac{\delta\,s} {\delta\,x}+c\frac{\delta\,s} {\delta\,x}}\right ]_{t+\Delta\,t}</math>
-where we follow the suggestion of Hindmarsh () and set <math>\alpha=1.5</math>.  The calculation is best explained in terms of the matrix equation '''Ax=b''', where the vector '''''x''''' (in vector '''answ''') contains the unknown thickness at the new time level; the vector '''''b''''' (in vector rhsd) contains the known thickness at the old time level, the known mass balance and the known flow at the old time level; and the matrix '''''A''''' contains coefficients relating to the thickness and surface elevation at the new time level (i.e., the finite difference approximations to the two terms within the bracketed quantities on the far right hand side).
-The algorithm proceeds in five stages.
-#Determine the diffusivities (D) from available thickness and surface elevation at the old time level; we approximate diffusivity at the new time level by the value at the old level.  Diffusivity is calculated on a staggered grid halfway between adjacent thickness points.  This is done by calling function '''[[getdiff]]''' and checking to see if the particular point under consideration is affected by a branch.
-#Determine the flow convergence (c) from available width and diffusivity at the old time level; again we approximate convergence at the new time level by the value at the old level. Diffusivity is calculated at the original thickness grid points.  This is done by calling function '''[[getconv]]''' and checking to see if the particular point under consideration is affected by a branch.
-#We are now in a position to fill matrix '''''A'''''.  We enter a value into the sparse form of A if the mask variable shows that this is appropriate (see '''function calcmask''').  This is done by repeated calls to '''putpcgc'''. In the case of a non-branching system, this results in a simple tridiagonal matrix.  A great deal of care is required to enter this information correctly for branched systems.  The known vector '''''b''''' is filled at the same time using the diffusive and convergence calculated above.
-#With the matrix '''''A''''' and vector '''''b''''' filled, we can obtain '''''x''''' by standard matrix inversion algorithms (in this case SLAP, see '''[[slapsolve]]'''). We only do this if there is a significant length of glacier etc.
-#Finally, unpack the new surface elevations from the matrix inversion output to the appropriate positions in the surface elevation array (usrf) using information in mask, and then find new thicknesses using bedrock elevation (topg) and glacier length.
+=== 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]]'''.<br>
-function getconv (contained in calcthck)
-This function returns the effect of flow convergence (i.e., varying widths) along the glacier.  This is determined from
-where standard symbols are defined in Table 1.  The term is scaled and is determined on the non-staggered grid that is used for thickness calculations.  It uses diffusivity from the staggered grid and widths on the non-staggered grid.  The function is called for each grid point in turn; this means that it is the responsibility of the calling routine to deal with branching in the system.
-This function is called from calcthck and the returned information is held in array conver.
-subroutine putpcgc (contained in calcthck)
-This simple subroutine fills a sparse matrix using the coefficients supplied.  The sparse matrix has a three column format (row, column and value) and can be used with the SLAP matrix algebra library.  The subroutine checks to ensure no zero values are entered.  The sparse matrix itself is held in module pcgdwk.
-This subroutine is called from calcthck and is called once per matrix element.
-subroutine slapsolv (contained in calcthck)
-This subroutine employs the SLAP matrix algebra library to solve the matrix inversion problem Ax=b, where A and b are known and x is required (in this case, ice thickness).  It expects a sparse form of the matrix to exist in module pcgdwk, and will first convert the anticipated three-column format into SLAP’s own storage format (ds2y).  The inversion itself is performed using one of a number of preconditioned conjugate-gradient methods available in SLAP (dslucs).  Finally, the convergence of the inversion is checked and the program stopped if convergence has not been achieved.
-This subroutine is called from calcthck.
-function calcflux
-This function returns the total horizontal flux of ice along the glacier.  This is determined from
-where standard symbols are defined in Table 1.  The flux is scaled and is determined on a staggered grid that has points positioned halfway between adjacent thickness grid points.  The algorithm takes into account any branching in the system.  If both up and downstream points have no ice then flux is set to zero.
-Note that this calculation is not needed for the main thickness-evolution calculation and is performed solely for diagnostic reasons. This function is called from valley and the returned information is held in array flux.
-function calcwdth
-This function returns the scaled width of the glacier.  This is determined using
-where standard symbols are defined in Table 1, and wr and Hr are (scaled) reference widths and thickness.  The parameter p is typically 0.5 and is hard coded in the function, while w1 is a small number incorporated to stop width becoming zero in the vent of zero ice thickness (which is important in the ice thickness solver).
-This function is called from valley and the returned information is held in array wdth.
-Module mass_balance (mass_balance.f90)
-function calcacab
-This function controls the calculation of the glacier’s surface mass balance.  It uses unscaled variables to do this except where variables have been scalled elsewhere (e.g., surface elevation).  t offer two mthods of determining ablation (day degree and simple energy balance) and has simple model for melt water retention as superimposed ice.
 The function is called once per year from the main valley model and has a daily time step.
-.	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.
+#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.
-.	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).
+#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).
-.	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).
+#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).
-.	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.
+#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.
-.	An estimation of the amount of precipitation falling as snow is then made based on diurnal air temperatures.
+#An estimation of the amount of precipitation falling as snow is then made based on diurnal air temperatures.
-.	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.
+#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.
-.	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.
+#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.
-.	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.
+#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 massbalseasons
-This function …
-function rainorsnw
-This function determines whether precipitation falls as rain or snow (it returns amount of rainfall) using air temperature.  If it is called with the flag not equal to 1, a simple comparison between daily mean temperature and a prescribed rain-snow transition temperature is made: if the daily mean is above the threshold then all precipitation falls as rain, otherwise none.
-If the function is called with flag set to 1, it estimates the rainfall by assuming a cosine for the diurnal variation of temperature about its mean ( ) with prescribed range ( from mean to peak NOT peak to peak)
-which can be rearranged to give time of day   (in radians) at which the threshold  is equaled
- .
-The proportion of the day for which temperatures exceed the threshold is therefore  .  If we assume that there is an equal probability of the precipitation falling at any time during the day then we simply need to multiply by the daily precipitation total to get the amount of rainfall.
-function finddegdays
-This function finds the number of positive day degrees ( ) during a day using air temperature.  If it is called with the flag not equal to 1, the number of day degrees is simply the daily mean temperature if its is above zero otherwise nothing.
-If the function is called with flag set to 1, it estimates the number of positive day degrees by assuming a cosine for the diurnal variation of temperature about its mean ( ) with prescribed range ( from mean to peak NOT peak to peak)
-which can be rearranged to give time of day   (in radians) at which the temperature equals zero
- .
-We now need to integrate   from the this time to midday ( ) and scale from radians into fractions of a day
-function lapsecalc
-This function simply applies the prescribed lapse rate to what ever variable is given.  This is usually done for air temperature (daily means but not amplitudes) and precipitation.  In the latter case, there is a guard in the calling routine that avoids creating precipitation if there was none in the base met. data.   The lapse rate is applied relative to the prescribed elevation of the met. data and everything is done in scaled units.
-function surfacemodel
-This function takes the estimated daily precipitation (total), rainfall and meltwater depths and evolves the accumulated snow, super-imposed ice and glacier ice stores.  We call firn the combined depths of snow and super-imposed ice.  The logic of the process is as follows.
-.	The depth of snow is increased by snowfall (total precipitation minus rainfall).
-.	The depth of super-imposed ice is increased by rainfall (assumed to refreeze).
-.	Use the total potential melt to melt snow.  If there is snow left then set melt to zero, otherwise reduce it by amount of snow melted.  Any melted snow is assumed to refreeze and become super-imposed ice.
-.	Estimate the depth of super-imposed ice that would result in the firn reaching its capacity.  Here firn is assumed to be the sum of the snow (s) and super-imposed ice (i) depths.  We assume that only a fixed proportion (w or wmax) of the firn can be super-imposed ice and estimate this depth as
-.	Any super-imposed ice above this capacity is assumed to have run off – or more accurately the melt water that would have formed super-imposed ice is rejected before it gets a chance to freeze.
-.	If any melt potential remains then proceed to melt super-imposed ice (reducing potential and increasing runoff accordingly).  If the melt potential exceeds the amount of super-imposed ice, remove all of and amend runoff and potential melt figures accordingly. Note that in the case of the degree day model, this requires a re-determination of the appropriate degree-day factor from that of snow to ice.
-.	If any melt potential remains then proceed to melt glacier ice ice (reducing potential and increasing runoff accordingly).
-Module energy_balance (energy_balance.f90)
-subroutine initialstart_ebm
-Sets up a number of parameters that are used in the energy balance calculation of surface melt.  Points to note include: scaling of time from second in day to radians.
-subroutine initial4day_ebm
-Calculate energy-budget parameters that vary daily.  In this case, this is the solar declination ( ) which is then used later on with the latitude ( ) to calculate the solar zenith (z) and azimuth ( )
-where   is the julian day expressed in radians, and
+=== Function calcmask ===
-the primes indicates that these are parameters used in the final calculation of zenith and azimuths.
+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:
-subroutine nexttimestep_ebm
+:*the point has ice thickness greater than zero;
-This is the main solver for the energy-balance model and steps both the surface temperature (T) and the accumulated melt forward in time.  The equation being solved is
+:*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.
-where   is the surface’s aspect and   its slope (in radians);   is the fraction of the day expressed in radians and relative to midday; A, B and C are free parameters;  is the co-albedo of the surface; d is the depth of the thermally-active layer,  its density and c its specific heat capacity (we currently assume snow throughout); the 2-m air temperature (T2m) is based on the daily mean and amplitude value supplied; and temperatures are in K.  Appropriate sets are taken to ensure that the sun does not sun when it is below the horizon.
+<br>
-We use an implicit time-marching scheme
-In these subroutines, time through the day is in radians.  If the new temperature is above the melt point then we estimate an amount of melt by converting energy with active layer into an equivalent melt depth.  We use the trapezium rule to do this.  There are two cases.  The first is where the preceeding temperature is at the melt point also; this is easy.  The second in when the preceeding temperature is below the melt point, in which case we need to estimate the fraction of the time step during which temperature was above the melt point (using linear interpolation).
+<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]]).''' </font>