<?xml version="1.0"?>
<feed xmlns="http://www.w3.org/2005/Atom" xml:lang="en">
	<id>https://source.geography.bristol.ac.uk/mediawiki/api.php?action=feedcontributions&amp;feedformat=atom&amp;user=RosDeath</id>
	<title>SourceWiki - User contributions [en]</title>
	<link rel="self" type="application/atom+xml" href="https://source.geography.bristol.ac.uk/mediawiki/api.php?action=feedcontributions&amp;feedformat=atom&amp;user=RosDeath"/>
	<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/wiki/Special:Contributions/RosDeath"/>
	<updated>2026-04-04T16:05:50Z</updated>
	<subtitle>User contributions</subtitle>
	<generator>MediaWiki 1.35.8</generator>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Model_mass_balance&amp;diff=4439</id>
		<title>Model mass balance</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Model_mass_balance&amp;diff=4439"/>
		<updated>2007-11-20T15:41:51Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: /* subroutine nexttimestep_ebm */&lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;There are two methods to calculate the mass balance within the model; &lt;br /&gt;
:*degree day method&lt;br /&gt;
:*simple energy balance model.&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
== Degree day method ==&lt;br /&gt;
The degree day method is the simplest method to calculate the ablation. In this type of model the potential ablation is determined as a linear function of the excess of daily mean air temperature above 0 °C. This potential ablation is transformed into an amount of melt through the application of degree day factors. The positive degree-day factor parameterises the physical processes that characterise the exchange of energy between the atmosphere and the surface. As a parameterisation of complex surface interchanges, degree-day factors are spatially and temporally variable. The degree day factor for snow and ice are two of the free parameters that are specified in '''[[Readoptimfile]]'''.&lt;br /&gt;
&lt;br /&gt;
=== function rainorsnw ===&lt;br /&gt;
&lt;br /&gt;
This function determines whether precipitation falls as rain or snow using air temperature(it returns amount of rainfall). If it is called with the flag not equal to 1, a simple comparison between the daily mean temperature and a prescribed rain-snow transition temperature is made: if the daily mean air temperature is above the threshold then all precipitation on that day falls as rain, otherwise it is assumed that all the precipitation is snow, so the function returns zero.&amp;lt;br&amp;gt;  &lt;br /&gt;
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 (&amp;lt;math&amp;gt;\bar{T}&amp;lt;/math&amp;gt;) with prescribed range, &amp;lt;math&amp;gt;\Delta\,T&amp;lt;/math&amp;gt;, (the range is from the mean to the peak, NOT peak to peak),&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;T(t)=\bar{T}-\Delta\,T \cos(t)&amp;lt;/math&amp;gt;&lt;br /&gt;
	 &lt;br /&gt;
which can be rearranged to give time of day &amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt; (in radians) at which the threshold &amp;lt;math&amp;gt;T^{'}&amp;lt;/math&amp;gt;  is reached, &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;t^{'}=\cos^{-1}\left({\frac{\bar{T}-T^{'}}{\Delta\,T}}\right)&amp;lt;/math&amp;gt; .&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
The proportion of the day for which temperatures exceed the threshold is therefore &amp;lt;math&amp;gt;(\pi\,-t^{'})/\pi\, &amp;lt;/math&amp;gt; .  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.&lt;br /&gt;
&lt;br /&gt;
=== function finddegdays ===&lt;br /&gt;
&lt;br /&gt;
This function finds the number of positive day degrees (''d'') 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, or none if the temperature is below zero.  &lt;br /&gt;
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 (&amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt;) with prescribed range, &amp;lt;math&amp;gt;\Delta\,T&amp;lt;/math&amp;gt;, (the range is from the mean to the peak, NOT peak to peak),&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;T(t)=\bar{T}-\Delta\,T \cos(t)&amp;lt;/math&amp;gt; &lt;br /&gt;
&lt;br /&gt;
which can be rearranged to give time of day, &amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt;, (in radians) at which the temperature equals zero, &amp;lt;br&amp;gt; &lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;t^{'}=\cos^{-1}\left({\frac{\bar{T}}{\Delta\,T}}\right)&amp;lt;/math&amp;gt; .&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
We now need to integrate ''T(t)'' from the this time to midday &amp;lt;math&amp;gt;\left({\pi\,}\right)&amp;lt;/math&amp;gt; and scale from radians into fractions of a day,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;d=\frac{1}{\pi\,} \int\limits_{t^{'}}^{\pi\,}\bar{T}-\Delta\,T\cos(t)dt&amp;lt;/math&amp;gt;&lt;br /&gt;
&lt;br /&gt;
:&amp;lt;math&amp;gt;=\frac{1}{\pi\,}\left[{\bar{T}\left({\pi\,-t^{'}}\right)+\Delta\,T\sin\left({t^{'}}\right)}\right]&amp;lt;/math&amp;gt;.&lt;br /&gt;
&lt;br /&gt;
=== function lapsecalc ===&lt;br /&gt;
&lt;br /&gt;
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 meteorological data.  &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
Although the free atmosphere temperature lapse rate is assumed to be constant at between – 6 °C to – 7 °C and this value is often used in studies investigating glacier mass balance response to climate change, the actual vertical gradient in surface temperatures is found to be spatially and temporally variable. The lapse rates are therefore opitmised in the model and so are defined in '''[[Readoptimfile]]''' . The lapse rate is applied relative to the elevation of the meteorological data, either the elevation of the meteorological station, or sea-level for ERA 40 data. All calculations are done in scaled units.&lt;br /&gt;
&lt;br /&gt;
=== function surfacemodel ===&lt;br /&gt;
&lt;br /&gt;
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.&lt;br /&gt;
# The depth of snow is increased by snowfall (total precipitation minus rainfall).&lt;br /&gt;
# The depth of super-imposed ice is increased by rainfall (assumed to refreeze).&lt;br /&gt;
# 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. &lt;br /&gt;
# 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,&lt;br /&gt;
&lt;br /&gt;
:::&amp;lt;math&amp;gt;i=w\left({s+i}\right)&amp;lt;/math&amp;gt; &lt;br /&gt;
&lt;br /&gt;
:::&amp;lt;math&amp;gt;i=\frac{ws}{1-w}&amp;lt;/math&amp;gt;&lt;br /&gt;
&amp;lt;ol&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;5&amp;quot;&amp;gt;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.&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;6&amp;quot;&amp;gt;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.&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;7&amp;quot;&amp;gt;If any melt potential remains then proceed to melt glacier ice ice (reducing potential and increasing runoff accordingly).&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;/ol&amp;gt;&lt;br /&gt;
&lt;br /&gt;
== Energy balance model ==&lt;br /&gt;
&lt;br /&gt;
=== subroutine initialstart_ebm ===&lt;br /&gt;
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.&lt;br /&gt;
&lt;br /&gt;
=== subroutine initial4day_ebm ===&lt;br /&gt;
&lt;br /&gt;
Calculates the energy-budget parameters that vary daily.  In this case, this is the solar declination (&amp;lt;math&amp;gt;\delta\,&amp;lt;/math&amp;gt;) which is then used later on with the latitude (&amp;lt;math&amp;gt;\phi\,&amp;lt;/math&amp;gt;) to calculate the solar zenith (z) and azimuth (&amp;lt;math&amp;gt;\psi\,&amp;lt;/math&amp;gt;), &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;&lt;br /&gt;
\begin{cases}&lt;br /&gt;
\delta\,= 0.006918 - 0.399912\cos\theta\,_{yr}+ 0.70257\sin\theta_{yr} -\\ 0.006758\cos2\theta_{yr} +0.000907\sin2\theta_{yr}- 0.002697\cos3\theta_{yr}+\\ 0.00148\sin3\theta_{yr} &lt;br /&gt;
\end{cases}&lt;br /&gt;
&amp;lt;/math&amp;gt;	 &lt;br /&gt;
&lt;br /&gt;
where &amp;lt;math&amp;gt;\theta_{yr}&amp;lt;/math&amp;gt; is the julian day expressed in radians. From this,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;z^{'}=&lt;br /&gt;
\begin{pmatrix}&lt;br /&gt;
  \sin\lambda\, &amp;amp; \sin\delta\, \\&lt;br /&gt;
  \cos\lambda\, &amp;amp; \cos\delta\,&lt;br /&gt;
\end{pmatrix}&lt;br /&gt;
&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
and,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;\psi^{'}=&lt;br /&gt;
\begin{pmatrix}&lt;br /&gt;
  \sin\lambda\, &amp;amp; \cos\delta\, \\&lt;br /&gt;
  \cos\lambda\, &amp;amp; \sin\delta\,&lt;br /&gt;
\end{pmatrix}&lt;br /&gt;
&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
the primes indicates that these are parameters used in the final calculation of zenith and azimuths.&lt;br /&gt;
&lt;br /&gt;
=== subroutine nexttimestep_ebm ===&lt;br /&gt;
&lt;br /&gt;
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, &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;\rho\,cd\frac{dT}{dt}=C\alpha\,^{'}I+BT_{2m}+AT &amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
::&amp;lt;math&amp;gt;I =I_{0}\left({\cos\,z\cos\phi\,+\sin\,z\sin\,\phi\,\cos\,\left({\psi\,-\beta\,}\right)}\right)&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
:&amp;lt;math&amp;gt;cos\,z=z^{'}_{1}+z^{'}_{2}\cos\,\theta_{d}&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
:&amp;lt;math&amp;gt;cos\,\psi\,=\frac{\psi\,^{'}_{1}\cos\,\theta\,_{d}-\psi\,^{'}_{2}}{\cos\,z}&amp;lt;/math&amp;gt;&lt;br /&gt;
&lt;br /&gt;
	 &lt;br /&gt;
&lt;br /&gt;
where &amp;lt;math&amp;gt;\rho\,&amp;lt;/math&amp;gt; is density, d is the depth of the thermally-active layer and c is the specific heat capacity (we currently assume snow throughout).  ''A'', ''B'' and ''C'' are free parameters. &amp;lt;math&amp;gt;\alpha\,^{'}&amp;lt;/math&amp;gt; is the co-albedo of the surface, ''I'' is the insolation at the surface; the 2-m air temperature (&amp;lt;math&amp;gt;T_{2m}&amp;lt;/math&amp;gt;) is based on the daily mean and amplitude value supplied. All temperatures are in K.&amp;lt;br&amp;gt;&lt;br /&gt;
&amp;lt;math&amp;gt;\beta\,&amp;lt;/math&amp;gt; is the surface’s aspect, &amp;lt;math&amp;gt;\phi\,&amp;lt;/math&amp;gt; is the slope (in radians)and &amp;lt;math&amp;gt;I_{0}&amp;lt;/math&amp;gt; is the insolation received at the top of the atmosphere, and is given a constant value of 1371 Wm&amp;lt;sup&amp;gt;-2&amp;lt;/sup&amp;gt;. &amp;lt;math&amp;gt;\theta_{d}\,&amp;lt;/math&amp;gt; is the fraction of the day expressed in radians and relative to midday. Appropriate sets are taken to ensure that there is no insolation when the sun is below the horizon. &amp;lt;br&amp;gt;&lt;br /&gt;
 &lt;br /&gt;
We use an implicit time-marching scheme, &lt;br /&gt;
&amp;lt;math&amp;gt;T^{n}-T^{n-1}=\frac{\delta\,t}{\rho\,cd}\left({C\alpha\,^{'}I^{n}+BT_{2m}^{n}+\frac{A}{2}\left({T^{n}-T^{n-1}}\right)}\right)&amp;lt;/math&amp;gt;.&lt;br /&gt;
&lt;br /&gt;
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 in the active layer into an equivalent melt depth.  We use the trapezium rule to do this.  There are two cases.  &lt;br /&gt;
#The first is where the preceeding temperature is at the melt point also; this is easy.  &lt;br /&gt;
#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).&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Model_mass_balance&amp;diff=4438</id>
		<title>Model mass balance</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Model_mass_balance&amp;diff=4438"/>
		<updated>2007-11-20T14:57:38Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: /* subroutine nexttimestep_ebm */&lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;There are two methods to calculate the mass balance within the model; &lt;br /&gt;
:*degree day method&lt;br /&gt;
:*simple energy balance model.&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
== Degree day method ==&lt;br /&gt;
The degree day method is the simplest method to calculate the ablation. In this type of model the potential ablation is determined as a linear function of the excess of daily mean air temperature above 0 °C. This potential ablation is transformed into an amount of melt through the application of degree day factors. The positive degree-day factor parameterises the physical processes that characterise the exchange of energy between the atmosphere and the surface. As a parameterisation of complex surface interchanges, degree-day factors are spatially and temporally variable. The degree day factor for snow and ice are two of the free parameters that are specified in '''[[Readoptimfile]]'''.&lt;br /&gt;
&lt;br /&gt;
=== function rainorsnw ===&lt;br /&gt;
&lt;br /&gt;
This function determines whether precipitation falls as rain or snow using air temperature(it returns amount of rainfall). If it is called with the flag not equal to 1, a simple comparison between the daily mean temperature and a prescribed rain-snow transition temperature is made: if the daily mean air temperature is above the threshold then all precipitation on that day falls as rain, otherwise it is assumed that all the precipitation is snow, so the function returns zero.&amp;lt;br&amp;gt;  &lt;br /&gt;
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 (&amp;lt;math&amp;gt;\bar{T}&amp;lt;/math&amp;gt;) with prescribed range, &amp;lt;math&amp;gt;\Delta\,T&amp;lt;/math&amp;gt;, (the range is from the mean to the peak, NOT peak to peak),&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;T(t)=\bar{T}-\Delta\,T \cos(t)&amp;lt;/math&amp;gt;&lt;br /&gt;
	 &lt;br /&gt;
which can be rearranged to give time of day &amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt; (in radians) at which the threshold &amp;lt;math&amp;gt;T^{'}&amp;lt;/math&amp;gt;  is reached, &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;t^{'}=\cos^{-1}\left({\frac{\bar{T}-T^{'}}{\Delta\,T}}\right)&amp;lt;/math&amp;gt; .&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
The proportion of the day for which temperatures exceed the threshold is therefore &amp;lt;math&amp;gt;(\pi\,-t^{'})/\pi\, &amp;lt;/math&amp;gt; .  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.&lt;br /&gt;
&lt;br /&gt;
=== function finddegdays ===&lt;br /&gt;
&lt;br /&gt;
This function finds the number of positive day degrees (''d'') 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, or none if the temperature is below zero.  &lt;br /&gt;
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 (&amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt;) with prescribed range, &amp;lt;math&amp;gt;\Delta\,T&amp;lt;/math&amp;gt;, (the range is from the mean to the peak, NOT peak to peak),&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;T(t)=\bar{T}-\Delta\,T \cos(t)&amp;lt;/math&amp;gt; &lt;br /&gt;
&lt;br /&gt;
which can be rearranged to give time of day, &amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt;, (in radians) at which the temperature equals zero, &amp;lt;br&amp;gt; &lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;t^{'}=\cos^{-1}\left({\frac{\bar{T}}{\Delta\,T}}\right)&amp;lt;/math&amp;gt; .&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
We now need to integrate ''T(t)'' from the this time to midday &amp;lt;math&amp;gt;\left({\pi\,}\right)&amp;lt;/math&amp;gt; and scale from radians into fractions of a day,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;d=\frac{1}{\pi\,} \int\limits_{t^{'}}^{\pi\,}\bar{T}-\Delta\,T\cos(t)dt&amp;lt;/math&amp;gt;&lt;br /&gt;
&lt;br /&gt;
:&amp;lt;math&amp;gt;=\frac{1}{\pi\,}\left[{\bar{T}\left({\pi\,-t^{'}}\right)+\Delta\,T\sin\left({t^{'}}\right)}\right]&amp;lt;/math&amp;gt;.&lt;br /&gt;
&lt;br /&gt;
=== function lapsecalc ===&lt;br /&gt;
&lt;br /&gt;
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 meteorological data.  &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
Although the free atmosphere temperature lapse rate is assumed to be constant at between – 6 °C to – 7 °C and this value is often used in studies investigating glacier mass balance response to climate change, the actual vertical gradient in surface temperatures is found to be spatially and temporally variable. The lapse rates are therefore opitmised in the model and so are defined in '''[[Readoptimfile]]''' . The lapse rate is applied relative to the elevation of the meteorological data, either the elevation of the meteorological station, or sea-level for ERA 40 data. All calculations are done in scaled units.&lt;br /&gt;
&lt;br /&gt;
=== function surfacemodel ===&lt;br /&gt;
&lt;br /&gt;
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.&lt;br /&gt;
# The depth of snow is increased by snowfall (total precipitation minus rainfall).&lt;br /&gt;
# The depth of super-imposed ice is increased by rainfall (assumed to refreeze).&lt;br /&gt;
# 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. &lt;br /&gt;
# 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,&lt;br /&gt;
&lt;br /&gt;
:::&amp;lt;math&amp;gt;i=w\left({s+i}\right)&amp;lt;/math&amp;gt; &lt;br /&gt;
&lt;br /&gt;
:::&amp;lt;math&amp;gt;i=\frac{ws}{1-w}&amp;lt;/math&amp;gt;&lt;br /&gt;
&amp;lt;ol&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;5&amp;quot;&amp;gt;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.&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;6&amp;quot;&amp;gt;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.&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;7&amp;quot;&amp;gt;If any melt potential remains then proceed to melt glacier ice ice (reducing potential and increasing runoff accordingly).&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;/ol&amp;gt;&lt;br /&gt;
&lt;br /&gt;
== Energy balance model ==&lt;br /&gt;
&lt;br /&gt;
=== subroutine initialstart_ebm ===&lt;br /&gt;
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.&lt;br /&gt;
&lt;br /&gt;
=== subroutine initial4day_ebm ===&lt;br /&gt;
&lt;br /&gt;
Calculates the energy-budget parameters that vary daily.  In this case, this is the solar declination (&amp;lt;math&amp;gt;\delta\,&amp;lt;/math&amp;gt;) which is then used later on with the latitude (&amp;lt;math&amp;gt;\phi\,&amp;lt;/math&amp;gt;) to calculate the solar zenith (z) and azimuth (&amp;lt;math&amp;gt;\psi\,&amp;lt;/math&amp;gt;), &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;&lt;br /&gt;
\begin{cases}&lt;br /&gt;
\delta\,= 0.006918 - 0.399912\cos\theta\,_{yr}+ 0.70257\sin\theta_{yr} -\\ 0.006758\cos2\theta_{yr} +0.000907\sin2\theta_{yr}- 0.002697\cos3\theta_{yr}+\\ 0.00148\sin3\theta_{yr} &lt;br /&gt;
\end{cases}&lt;br /&gt;
&amp;lt;/math&amp;gt;	 &lt;br /&gt;
&lt;br /&gt;
where &amp;lt;math&amp;gt;\theta_{yr}&amp;lt;/math&amp;gt; is the julian day expressed in radians. From this,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;z^{'}=&lt;br /&gt;
\begin{pmatrix}&lt;br /&gt;
  \sin\lambda\, &amp;amp; \sin\delta\, \\&lt;br /&gt;
  \cos\lambda\, &amp;amp; \cos\delta\,&lt;br /&gt;
\end{pmatrix}&lt;br /&gt;
&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
and,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;\psi^{'}=&lt;br /&gt;
\begin{pmatrix}&lt;br /&gt;
  \sin\lambda\, &amp;amp; \cos\delta\, \\&lt;br /&gt;
  \cos\lambda\, &amp;amp; \sin\delta\,&lt;br /&gt;
\end{pmatrix}&lt;br /&gt;
&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
the primes indicates that these are parameters used in the final calculation of zenith and azimuths.&lt;br /&gt;
&lt;br /&gt;
=== subroutine nexttimestep_ebm ===&lt;br /&gt;
&lt;br /&gt;
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, &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;\rho\,cd\frac{dT}{dt}=C\alpha\,^{'}I+BT_{2m}+AT &amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
::&amp;lt;math&amp;gt;I =I_{0}\left({\cos\,z\cos\phi\,+\sin\,z\sin\,\phi\,\cos\,\left({\psi\,-\beta\,}\right)}\right)&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
:&amp;lt;math&amp;gt;cos\,z=z^{'}_{1}+z^{'}_{2}\cos\,\theta_{d}&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
:&amp;lt;math&amp;gt;cos\,\psi\,=\frac{\psi\,^{'}_{1}\cos\,\theta\,_{d}-\psi\,^{'}_{2}}{\cos\,z}&amp;lt;/math&amp;gt;&lt;br /&gt;
&lt;br /&gt;
	 &lt;br /&gt;
&lt;br /&gt;
where &amp;lt;math&amp;gt;\rho\,&amp;lt;/math&amp;gt; is density, d is the depth of the thermally-active layer and c is the specific heat capacity (we currently assume snow throughout).  ''A'', ''B'' and ''C'' are free parameters. &amp;lt;math&amp;gt;\alpha\,^{'}&amp;lt;/math&amp;gt; is the co-albedo of the surface, ''I'' is the insolation at the surface; the 2-m air temperature (&amp;lt;math&amp;gt;T_{2m}&amp;lt;/math&amp;gt;) is based on the daily mean and amplitude value supplied. All temperatures are in K.&amp;lt;br&amp;gt;&lt;br /&gt;
&amp;lt;math&amp;gt;\beta\,&amp;lt;/math&amp;gt; is the surface’s aspect, &amp;lt;math&amp;gt;\phi\,&amp;lt;/math&amp;gt; is the slope (in radians)and &amp;lt;math&amp;gt;I_{0}&amp;lt;/math&amp;gt; is the insolation received at the top of the atmosphere, and is given a constant value of 1371 Wm&amp;lt;sup&amp;gt;2&amp;lt;/sup&amp;gt;. &amp;lt;math&amp;gt;\theta_{d}\,&amp;lt;/math&amp;gt; is the fraction of the day expressed in radians and relative to midday. Appropriate sets are taken to ensure that there is no insolation when the sun is below the horizon. &amp;lt;br&amp;gt;&lt;br /&gt;
 &lt;br /&gt;
We use an implicit time-marching scheme, &lt;br /&gt;
&amp;lt;math&amp;gt;T^{n}-T^{n-1}=\frac{\delta\,t}{\rho\,cd}\left({C\alpha\,^{'}I^{n}+BT_{2m}^{n}+\frac{A}{2}\left({T^{n}-T^{n-1}}\right)}\right)&amp;lt;/math&amp;gt;.&lt;br /&gt;
&lt;br /&gt;
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 in the active layer into an equivalent melt depth.  We use the trapezium rule to do this.  There are two cases.  &lt;br /&gt;
#The first is where the preceeding temperature is at the melt point also; this is easy.  &lt;br /&gt;
#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).&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4437</id>
		<title>Phase space mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4437"/>
		<updated>2007-11-20T14:20:39Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: &lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;In this mode the model explores the phase space of the various parameters that are being optimised. The parameter values set in '''[[Readoptimfile]]''' have boundary values within which parameter values may vary. The subroutine '''Space_survey''' runs the model over the range of possible parameter values. For each combination of parameter values a cost function is determined,&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;C = \frac{1}{n_{i}}\sum_{i} \left|{M_{i}-m_{i}}\right|+\frac{1}{n_{j}}\sum_{j} \left|{M_{j}-m_{j}}\right|+\frac{1}{n_{k}}\sum_{k} \left|{M_{k}-m_{k}}\right|+...&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
where, C is the cost function value, M is the measured value, m is the model predicted value, n is the normalisation constant, and the suffix i, j, k, refers to the different optimisation datasets.&lt;br /&gt;
&lt;br /&gt;
The output from this mode of model run shows the surface of the phase space. In the example below, the response surface is shown for the parameters air temperature lapse rate (&amp;amp;deg;C km&amp;lt;sup&amp;gt;-1&amp;lt;/sup&amp;gt;) and precipitation lapse rate (m m&amp;lt;sup&amp;gt;-1&amp;lt;/sup&amp;gt;day&amp;lt;sup&amp;gt;-1&amp;lt;/sup&amp;gt;). The colour of the surface corresponds to the cost function for each combination of temperature and precipitation lapse rates.&amp;lt;br&amp;gt;&lt;br /&gt;
[[Image:phasespace_ex.jpg]]&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4436</id>
		<title>Phase space mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4436"/>
		<updated>2007-11-20T14:20:14Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: &lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;In this mode the model explores the phase space of the various parameters that are being optimised. The parameter values set in '''[[Readoptimfile]]''' have boundary values within which parameter values may vary. The subroutine '''Space_survey''' runs the model over the range of possible parameter values. For each combination of parameter values a cost function is determined,&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;C = \frac{1}{n_{i}}\sum_{i} \left|{M_{i}-m_{i}}\right|+\frac{1}{n_{j}}\sum_{j} \left|{M_{j}-m_{j}}\right|+\frac{1}{n_{k}}\sum_{k} \left|{M_{k}-m_{k}}\right|+...&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
where, C is the cost function value, M is the measured value, m is the model predicted value, n is the normalisation constant, and the suffix i, j, k, refers to the different optimisation datasets.&lt;br /&gt;
&lt;br /&gt;
The output from this mode of model run shows the surface of the phase space. In the example below, the response surface is shown for the parameters air temperature lapse rate (&amp;amp;deg;C km&amp;lt;sup&amp;gt;-1&amp;lt;/sup&amp;gt;) and precipitation lapse rate (m m&amp;lt;sup&amp;gt;-1&amp;lt;/sup&amp;gt;day&amp;lt;sup&amp;gt;-1&amp;lt;/sup&amp;gt;&amp;lt;br&amp;gt;. The colour of the surface corresponds to the cost function for each combination of temperature and precipitation lapse rates.&lt;br /&gt;
[[Image:phasespace_ex.jpg]]&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4435</id>
		<title>Phase space mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4435"/>
		<updated>2007-11-20T14:18:33Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: &lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;In this mode the model explores the phase space of the various parameters that are being optimised. The parameter values set in '''[[Readoptimfile]]''' have boundary values within which parameter values may vary. The subroutine '''Space_survey''' runs the model over the range of possible parameter values. For each combination of parameter values a cost function is determined,&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;C = \frac{1}{n_{i}}\sum_{i} \left|{M_{i}-m_{i}}\right|+\frac{1}{n_{j}}\sum_{j} \left|{M_{j}-m_{j}}\right|+\frac{1}{n_{k}}\sum_{k} \left|{M_{k}-m_{k}}\right|+...&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
where, C is the cost function value, M is the measured value, m is the model predicted value, n is the normalisation constant, and the suffix i, j, k, refers to the different optimisation datasets.&lt;br /&gt;
&lt;br /&gt;
The output from this mode of model run shows the surface of the phase space. In the example below, the response surface when optimising for the parameters air temperature lapse rate (&amp;amp;deg;C km&amp;lt;sup&amp;gt;-1&amp;lt;/sup&amp;gt;) and precipitation lapse rate (m m&amp;lt;sup&amp;gt;-1&amp;lt;/sup&amp;gt;day&amp;lt;sup&amp;gt;-1&amp;lt;/sup&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
[[Image:phasespace_ex.jpg]]&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4434</id>
		<title>Phase space mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4434"/>
		<updated>2007-11-20T12:01:53Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: &lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;In this mode the model explores the phase space of the various parameters that are being optimised. The parameter values set in '''[[Readoptimfile]]''' have boundary values within which parameter values may vary. The subroutine '''Space_survey''' runs the model over the range of possible parameter values. For each combination of parameter values a cost function is determined,&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;C = \frac{1}{n_{i}}\sum_{i} \left|{M_{i}-m_{i}}\right|+\frac{1}{n_{j}}\sum_{j} \left|{M_{j}-m_{j}}\right|+\frac{1}{n_{k}}\sum_{k} \left|{M_{k}-m_{k}}\right|+...&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
where, C is the cost function value, M is the measured value, m is the model predicted value, n is the normalisation constant, and the suffix i, j, k, refers to the different optimisation datasets.&lt;br /&gt;
&lt;br /&gt;
The output from this mode of model run shows the surface of the phase space. In the example below, the response surface when optimising for the parameters air temperature lapse rate (&amp;amp;deg;C km&amp;lt;sup&amp;gt;-1&amp;lt;/sup&amp;gt;) and precipitation lapse rate (mm&amp;lt;br&amp;gt;&lt;br /&gt;
[[Image:phasespace_ex.jpg]]&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4433</id>
		<title>Phase space mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4433"/>
		<updated>2007-11-19T18:18:44Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: &lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;In this mode the model explores the phase space of the various parameters that are being optimised. The parameter values set in '''[[Readoptimfile]]''' have boundary values within which parameter values may vary. The subroutine '''Space_survey''' runs the model over the range of possible parameter values. For each combination of parameter values a cost function is determined,&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;C = \frac{1}{n_{i}}\sum_{i} \left|{M_{i}-m_{i}}\right|+\frac{1}{n_{j}}\sum_{j} \left|{M_{j}-m_{j}}\right|+\frac{1}{n_{k}}\sum_{k} \left|{M_{k}-m_{k}}\right|+...&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
where, C is the cost function value, M is the measured value, m is the model predicted value, n is the normalisation constant, and the suffix i, j, k, refers to the different optimisation datasets.&lt;br /&gt;
&lt;br /&gt;
The output from this mode of model run shows the surface of the phase space. In the example below, the response surface when optimising for the parameters air temperature lapse rate (&amp;amp;deg;C km&amp;lt;sup&amp;gt;-1&amp;lt;/sup&amp;gt;) and precipitation lapse rate (mm[[Image:phasespace_ex.jpg]]&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4432</id>
		<title>Phase space mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4432"/>
		<updated>2007-11-19T18:07:06Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: &lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;In this mode the model explores the phase space of the various parameters that are being optimised. The parameter values set in '''[[Readoptimfile]]''' have boundary values within which parameter values may vary. The subroutine '''Space_survey''' runs the model over the range of possible parameter values. For each combination of parameter values a cost function is determined,&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;C = \frac{1}{n_{i}}\sum_{i} \left|{M_{i}-m_{i}}\right|+\frac{1}{n_{j}}\sum_{j} \left|{M_{j}-m_{j}}\right|+\frac{1}{n_{k}}\sum_{k} \left|{M_{k}-m_{k}}\right|+...&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
where, C is the cost function value, M is the measured value, m is the model predicted value, n is the normalisation constant, and the suffix i, j, k, refers to the different optimisation datasets.&lt;br /&gt;
&lt;br /&gt;
The output from this mode of model run shows the surface of the phase space. In the example below, the response surface for optimising the parameters air temperature lapse rate (&amp;amp;deg;C km&amp;lt;sup&amp;gt;-1&amp;lt;/sup&amp;gt;)[[Image:phasespace_ex.jpg]]&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4431</id>
		<title>Phase space mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4431"/>
		<updated>2007-11-19T18:06:57Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: &lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;In this mode the model explores the phase space of the various parameters that are being optimised. The parameter values set in '''[[Readoptimfile]]''' have boundary values within which parameter values may vary. The subroutine '''Space_survey''' runs the model over the range of possible parameter values. For each combination of parameter values a cost function is determined,&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;C = \frac{1}{n_{i}}\sum_{i} \left|{M_{i}-m_{i}}\right|+\frac{1}{n_{j}}\sum_{j} \left|{M_{j}-m_{j}}\right|+\frac{1}{n_{k}}\sum_{k} \left|{M_{k}-m_{k}}\right|+...&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
where, C is the cost function value, M is the measured value, m is the model predicted value, n is the normalisation constant, and the suffix i, j, k, refers to the different optimisation datasets.&lt;br /&gt;
&lt;br /&gt;
The output from this mode of model run shows the surface of the phase space. In the example below, the response surface for optimising the parameters air temperature lapse rate (&amp;amp;deg;C km&amp;lt;sup&amp;gt;-1&amp;lt;/sup)[[Image:phasespace_ex.jpg]]&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4430</id>
		<title>Phase space mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4430"/>
		<updated>2007-11-19T18:05:54Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: &lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;In this mode the model explores the phase space of the various parameters that are being optimised. The parameter values set in '''[[Readoptimfile]]''' have boundary values within which parameter values may vary. The subroutine '''Space_survey''' runs the model over the range of possible parameter values. For each combination of parameter values a cost function is determined,&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;C = \frac{1}{n_{i}}\sum_{i} \left|{M_{i}-m_{i}}\right|+\frac{1}{n_{j}}\sum_{j} \left|{M_{j}-m_{j}}\right|+\frac{1}{n_{k}}\sum_{k} \left|{M_{k}-m_{k}}\right|+...&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
where, C is the cost function value, M is the measured value, m is the model predicted value, n is the normalisation constant, and the suffix i, j, k, refers to the different optimisation datasets.&lt;br /&gt;
&lt;br /&gt;
The output from this mode of model run shows the surface of the phase space. In the example below, the response surface for optimising the parameters air temperature lapse rate (&amp;amp;deg;C km&amp;lt;sup&amp;gt;-1)[[Image:phasespace_ex.jpg]]&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4429</id>
		<title>Phase space mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4429"/>
		<updated>2007-11-19T18:05:36Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: &lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;In this mode the model explores the phase space of the various parameters that are being optimised. The parameter values set in '''[[Readoptimfile]]''' have boundary values within which parameter values may vary. The subroutine '''Space_survey''' runs the model over the range of possible parameter values. For each combination of parameter values a cost function is determined,&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;C = \frac{1}{n_{i}}\sum_{i} \left|{M_{i}-m_{i}}\right|+\frac{1}{n_{j}}\sum_{j} \left|{M_{j}-m_{j}}\right|+\frac{1}{n_{k}}\sum_{k} \left|{M_{k}-m_{k}}\right|+...&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
where, C is the cost function value, M is the measured value, m is the model predicted value, n is the normalisation constant, and the suffix i, j, k, refers to the different optimisation datasets.&lt;br /&gt;
&lt;br /&gt;
The output from this mode of model run shows the surface of the phase space. In the example below, the response surface for optimising the parameters air temperature lapse rate (&amp;amp;deg;C km&amp;lt;sup&amp;gt;-1 )[[Image:phasespace_ex.jpg]]&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4428</id>
		<title>Phase space mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4428"/>
		<updated>2007-11-19T18:05:26Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: &lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;In this mode the model explores the phase space of the various parameters that are being optimised. The parameter values set in '''[[Readoptimfile]]''' have boundary values within which parameter values may vary. The subroutine '''Space_survey''' runs the model over the range of possible parameter values. For each combination of parameter values a cost function is determined,&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;C = \frac{1}{n_{i}}\sum_{i} \left|{M_{i}-m_{i}}\right|+\frac{1}{n_{j}}\sum_{j} \left|{M_{j}-m_{j}}\right|+\frac{1}{n_{k}}\sum_{k} \left|{M_{k}-m_{k}}\right|+...&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
where, C is the cost function value, M is the measured value, m is the model predicted value, n is the normalisation constant, and the suffix i, j, k, refers to the different optimisation datasets.&lt;br /&gt;
&lt;br /&gt;
The output from this mode of model run shows the surface of the phase space. In the example below, the response surface for optimising the parameters air temperature lapse rate (&amp;amp;deg;C km&amp;lt;sup&amp;gt;-1)[[Image:phasespace_ex.jpg]]&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4427</id>
		<title>Phase space mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4427"/>
		<updated>2007-11-19T18:05:04Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: &lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;In this mode the model explores the phase space of the various parameters that are being optimised. The parameter values set in '''[[Readoptimfile]]''' have boundary values within which parameter values may vary. The subroutine '''Space_survey''' runs the model over the range of possible parameter values. For each combination of parameter values a cost function is determined,&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;C = \frac{1}{n_{i}}\sum_{i} \left|{M_{i}-m_{i}}\right|+\frac{1}{n_{j}}\sum_{j} \left|{M_{j}-m_{j}}\right|+\frac{1}{n_{k}}\sum_{k} \left|{M_{k}-m_{k}}\right|+...&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
where, C is the cost function value, M is the measured value, m is the model predicted value, n is the normalisation constant, and the suffix i, j, k, refers to the different optimisation datasets.&lt;br /&gt;
&lt;br /&gt;
The output from this mode of model run shows the surface of the phase space. In the example below, the response surface for optimising the parameters air temperature lapse rate (&amp;amp;deg;C km&amp;lt;sup&amp;gt;-1[[Image:phasespace_ex.jpg]]&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4426</id>
		<title>Phase space mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4426"/>
		<updated>2007-11-19T17:59:54Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: &lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;In this mode the model explores the phase space of the various parameters that are being optimised. The parameter values set in '''[[Readoptimfile]]''' have boundary values within which parameter values may vary. The subroutine '''Space_survey''' runs the model over the range of possible parameter values. For each combination of parameter values a cost function is determined,&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;C = \frac{1}{n_{i}}\sum_{i} \left|{M_{i}-m_{i}}\right|+\frac{1}{n_{j}}\sum_{j} \left|{M_{j}-m_{j}}\right|+\frac{1}{n_{k}}\sum_{k} \left|{M_{k}-m_{k}}\right|+...&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
where, C is the cost function value, M is the measured value, m is the model predicted value, n is the normalisation constant, and the suffix i, j, k, refers to the different optimisation datasets.&lt;br /&gt;
&lt;br /&gt;
The output from this mode of model run shows the surface of the phase space. In the example below, the response surface for optimising the parameters air temperature lapse rate (&amp;amp;deg;C km [[Image:phasespace_ex.jpg]]&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4425</id>
		<title>Phase space mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4425"/>
		<updated>2007-11-19T17:54:46Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: &lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;In this mode the model explores the phase space of the various parameters that are being optimised. The parameter values set in '''[[Readoptimfile]]''' have boundary values within which parameter values may vary. The subroutine '''Space_survey''' runs the model over the range of possible parameter values. For each combination of parameter values a cost function is determined,&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;C = \frac{1}{n_{i}}\sum_{i} \left|{M_{i}-m_{i}}\right|+\frac{1}{n_{j}}\sum_{j} \left|{M_{j}-m_{j}}\right|+\frac{1}{n_{k}}\sum_{k} \left|{M_{k}-m_{k}}\right|+...&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
where, C is the cost function value, M is the measured value, m is the model predicted value, n is the normalisation constant, and the suffix i, j, k, refers to the different optimisation datasets.&lt;br /&gt;
&lt;br /&gt;
The output from this mode of model run shows the surface of the phase space. In the example below, the response surface for optimising the parameters air temperature lapse rate (&amp;lt;math&amp;gt;1^{\circ} C km^{-1}&amp;lt;/math&amp;gt;) [[Image:phasespace_ex.jpg]](&amp;amp;deg;C km&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4424</id>
		<title>Phase space mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4424"/>
		<updated>2007-11-19T17:45:58Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: &lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;In this mode the model explores the phase space of the various parameters that are being optimised. The parameter values set in '''[[Readoptimfile]]''' have boundary values within which parameter values may vary. The subroutine '''Space_survey''' runs the model over the range of possible parameter values. For each combination of parameter values a cost function is determined,&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;C = \frac{1}{n_{i}}\sum_{i} \left|{M_{i}-m_{i}}\right|+\frac{1}{n_{j}}\sum_{j} \left|{M_{j}-m_{j}}\right|+\frac{1}{n_{k}}\sum_{k} \left|{M_{k}-m_{k}}\right|+...&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
where, C is the cost function value, M is the measured value, m is the model predicted value, n is the normalisation constant, and the suffix i, j, k, refers to the different optimisation datasets.&lt;br /&gt;
&lt;br /&gt;
The output from this mode of model run shows the surface of the phase space. In the example below, the response surface for optimising the parameters air temperature lapse rate (&amp;lt;math&amp;gt;1^{\circ} C km^{-1}&amp;lt;/math&amp;gt;) [[Image:phasespace_ex.jpg]]&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4423</id>
		<title>Phase space mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4423"/>
		<updated>2007-11-19T17:45:42Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: &lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;In this mode the model explores the phase space of the various parameters that are being optimised. The parameter values set in '''[[Readoptimfile]]''' have boundary values within which parameter values may vary. The subroutine '''Space_survey''' runs the model over the range of possible parameter values. For each combination of parameter values a cost function is determined,&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;C = \frac{1}{n_{i}}\sum_{i} \left|{M_{i}-m_{i}}\right|+\frac{1}{n_{j}}\sum_{j} \left|{M_{j}-m_{j}}\right|+\frac{1}{n_{k}}\sum_{k} \left|{M_{k}-m_{k}}\right|+...&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
where, C is the cost function value, M is the measured value, m is the model predicted value, n is the normalisation constant, and the suffix i, j, k, refers to the different optimisation datasets.&lt;br /&gt;
&lt;br /&gt;
The output from this mode of model run shows the surface of the phase space. In the example below, the response surface for optimising the parameters air temperature lapse rate (&amp;lt;math&amp;gt;^{\circ} C km^{-1}&amp;lt;/math&amp;gt;) [[Image:phasespace_ex.jpg]]&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4422</id>
		<title>Phase space mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4422"/>
		<updated>2007-11-19T17:45:28Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: &lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;In this mode the model explores the phase space of the various parameters that are being optimised. The parameter values set in '''[[Readoptimfile]]''' have boundary values within which parameter values may vary. The subroutine '''Space_survey''' runs the model over the range of possible parameter values. For each combination of parameter values a cost function is determined,&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;C = \frac{1}{n_{i}}\sum_{i} \left|{M_{i}-m_{i}}\right|+\frac{1}{n_{j}}\sum_{j} \left|{M_{j}-m_{j}}\right|+\frac{1}{n_{k}}\sum_{k} \left|{M_{k}-m_{k}}\right|+...&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
where, C is the cost function value, M is the measured value, m is the model predicted value, n is the normalisation constant, and the suffix i, j, k, refers to the different optimisation datasets.&lt;br /&gt;
&lt;br /&gt;
The output from this mode of model run shows the surface of the phase space. In the example below, the response surface for optimising the parameters air temperature lapse rate (&amp;lt;math&amp;gt;\^{circ} C km^{-1}&amp;lt;/math&amp;gt;) [[Image:phasespace_ex.jpg]]&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4421</id>
		<title>Phase space mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4421"/>
		<updated>2007-11-19T17:45:07Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: &lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;In this mode the model explores the phase space of the various parameters that are being optimised. The parameter values set in '''[[Readoptimfile]]''' have boundary values within which parameter values may vary. The subroutine '''Space_survey''' runs the model over the range of possible parameter values. For each combination of parameter values a cost function is determined,&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;C = \frac{1}{n_{i}}\sum_{i} \left|{M_{i}-m_{i}}\right|+\frac{1}{n_{j}}\sum_{j} \left|{M_{j}-m_{j}}\right|+\frac{1}{n_{k}}\sum_{k} \left|{M_{k}-m_{k}}\right|+...&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
where, C is the cost function value, M is the measured value, m is the model predicted value, n is the normalisation constant, and the suffix i, j, k, refers to the different optimisation datasets.&lt;br /&gt;
&lt;br /&gt;
The output from this mode of model run shows the surface of the phase space. In the example below, the response surface for optimising the parameters air temperature lapse rate (&amp;lt;math&amp;gt;\circ C km^{-1}&amp;lt;/math&amp;gt;) [[Image:phasespace_ex.jpg]]&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4420</id>
		<title>Phase space mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4420"/>
		<updated>2007-11-19T17:44:56Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: &lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;In this mode the model explores the phase space of the various parameters that are being optimised. The parameter values set in '''[[Readoptimfile]]''' have boundary values within which parameter values may vary. The subroutine '''Space_survey''' runs the model over the range of possible parameter values. For each combination of parameter values a cost function is determined,&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;C = \frac{1}{n_{i}}\sum_{i} \left|{M_{i}-m_{i}}\right|+\frac{1}{n_{j}}\sum_{j} \left|{M_{j}-m_{j}}\right|+\frac{1}{n_{k}}\sum_{k} \left|{M_{k}-m_{k}}\right|+...&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
where, C is the cost function value, M is the measured value, m is the model predicted value, n is the normalisation constant, and the suffix i, j, k, refers to the different optimisation datasets.&lt;br /&gt;
&lt;br /&gt;
The output from this mode of model run shows the surface of the phase space. In the example below, the response surface for optimising the parameters air temperature lapse rate (&amp;lt;math&amp;gt;\circC km^{-1}&amp;lt;/math&amp;gt;) [[Image:phasespace_ex.jpg]]&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4419</id>
		<title>Phase space mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4419"/>
		<updated>2007-11-19T17:32:52Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: &lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;In this mode the model explores the phase space of the various parameters that are being optimised. The parameter values set in '''[[Readoptimfile]]''' have boundary values within which parameter values may vary. The subroutine '''Space_survey''' runs the model over the range of possible parameter values. For each combination of parameter values a cost function is determined,&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;C = \frac{1}{n_{i}}\sum_{i} \left|{M_{i}-m_{i}}\right|+\frac{1}{n_{j}}\sum_{j} \left|{M_{j}-m_{j}}\right|+\frac{1}{n_{k}}\sum_{k} \left|{M_{k}-m_{k}}\right|+...&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
where, C is the cost function value, M is the measured value, m is the model predicted value, n is the normalisation constant, and the suffix i, j, k, refers to the different optimisation datasets.&lt;br /&gt;
&lt;br /&gt;
The output from this mode of model run shows the surface of the phase space. In the example below, the response surface for optimising the parameters air temperature lapse rate (&amp;lt;math&amp;gt;\degC km^{-1}&amp;lt;/math&amp;gt;) [[Image:phasespace_ex.jpg]]&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=File:Phasespace_ex.jpg&amp;diff=4418</id>
		<title>File:Phasespace ex.jpg</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=File:Phasespace_ex.jpg&amp;diff=4418"/>
		<updated>2007-11-19T17:19:40Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: &lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4417</id>
		<title>Phase space mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4417"/>
		<updated>2007-11-19T17:19:19Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: &lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;In this mode the model explores the phase space of the various parameters that are being optimised. The parameter values set in '''[[Readoptimfile]]''' have boundary values within which parameter values may vary. The subroutine '''Space_survey''' runs the model over the range of possible parameter values. For each combination of parameter values a cost function is determined,&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;C = \frac{1}{n_{i}}\sum_{i} \left|{M_{i}-m_{i}}\right|+\frac{1}{n_{j}}\sum_{j} \left|{M_{j}-m_{j}}\right|+\frac{1}{n_{k}}\sum_{k} \left|{M_{k}-m_{k}}\right|+...&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
where, C is the cost function value, M is the measured value, m is the model predicted value, n is the normalisation constant, and the suffix i, j, k, refers to the different optimisation datasets.&lt;br /&gt;
&lt;br /&gt;
The output from this mode of model run shows the surface of the phase space.[[Image:phasespace_ex.jpg]]&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=File:Phasespace_ex1.JPG&amp;diff=4416</id>
		<title>File:Phasespace ex1.JPG</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=File:Phasespace_ex1.JPG&amp;diff=4416"/>
		<updated>2007-11-19T16:37:47Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: &lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;file that shows the phase space.&lt;br /&gt;
&lt;br /&gt;
back to [[Skadia]].&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=File:Phasespace_ex1.JPG&amp;diff=4415</id>
		<title>File:Phasespace ex1.JPG</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=File:Phasespace_ex1.JPG&amp;diff=4415"/>
		<updated>2007-11-19T16:22:12Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: &lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;file that shows the phase space&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=File:Phasespace_ex1.JPG&amp;diff=4414</id>
		<title>File:Phasespace ex1.JPG</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=File:Phasespace_ex1.JPG&amp;diff=4414"/>
		<updated>2007-11-19T16:21:13Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: file that shows the pahse space&lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;file that shows the pahse space&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4413</id>
		<title>Phase space mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4413"/>
		<updated>2007-11-16T18:32:51Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: &lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;In this mode the model explores the phase space of the various parameters that are being optimised. The parameter values set in '''[[Readoptimfile]]''' have boundary values within which parameter values may vary. The subroutine '''Space_survey''' runs the model over the range of possible parameter values. For each combination of parameter values a cost function is determined,&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;C = \frac{1}{n_{i}}\sum_{i} \left|{M_{i}-m_{i}}\right|+\frac{1}{n_{j}}\sum_{j} \left|{M_{j}-m_{j}}\right|+\frac{1}{n_{k}}\sum_{k} \left|{M_{k}-m_{k}}\right|+...&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
where, C is the cost function value, M is the measured value, m is the model predicted value, n is the normalisation constant, and the suffix i, j, k, refers to the different optimisation datasets.&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4412</id>
		<title>Phase space mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4412"/>
		<updated>2007-11-16T18:32:31Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: &lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;In this mode the model explores the phase space of the various parameters that are being optimised. The parameter values set in '''[[Readoptimfile]]''' have boundary values within which parameter values may vary. The subroutine Space_survey runs the model over the range of possible parameter values. For each combination of parameter values a cost function is determined,&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;C = \frac{1}{n_{i}}\sum_{i} \left|{M_{i}-m_{i}}\right|+\frac{1}{n_{j}}\sum_{j} \left|{M_{j}-m_{j}}\right|+\frac{1}{n_{k}}\sum_{k} \left|{M_{k}-m_{k}}\right|+...&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
where, C is the cost function value, M is the measured value, m is the model predicted value, n is the normalisation constant, and the suffix i, j, k, refers to the different optimisation datasets.&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4411</id>
		<title>Phase space mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4411"/>
		<updated>2007-11-16T18:31:57Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: &lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;In this mode the model explores the phase space of the various parameters that are being optimised. The parameter values set in '''[[Readoptimfile]]''' have boundary values within which parameter values may vary. The subroutine Space_survey runs the model over the range of possible parameter values. For each combination of parameter values a cost function is determined,&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4410</id>
		<title>Phase space mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Phase_space_mode&amp;diff=4410"/>
		<updated>2007-11-16T18:25:03Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: &lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;In this mode the model explores the phase space of the various parameters that are being optimised. The parameter values set in [[Readoptimfile]] have boundary values within which parameter values may vary. The subroutine Space_survey runs the model over the range of possible parameter values. For each combination of parameter values a cost function is determined,&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Optimisation_mode&amp;diff=4409</id>
		<title>Optimisation mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Optimisation_mode&amp;diff=4409"/>
		<updated>2007-11-16T17:48:38Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: /* Optimisation mode */&lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;== Optimisation mode ==&lt;br /&gt;
In this mode, the genetic algorithm solver is employed to determine the best values for the free parameters so that the difference between the modelled and the observed glacier response to climatic forcing is minimized. The cost function is determined by,&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;C = \frac{1}{n_{i}}\sum_{i} \left|{M_{i}-m_{i}}\right|+\frac{1}{n_{j}}\sum_{j} \left|{M_{j}-m_{j}}\right|+\frac{1}{n_{k}}\sum_{k} \left|{M_{k}-m_{k}}\right|+...&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
where, C is the cost function value, M is the measured value, m is the model predicted value, n is the normalisation constant, and the suffix i, j, k, refers to the different optimisation datasets. &lt;br /&gt;
&lt;br /&gt;
To run in this mode, the subroutine '''Creategadriver''' is called to create the Ga.inp file, which takes the parameters that are going to be optimised and puts them in a format that the genetic algorithm solver can read. &lt;br /&gt;
The genetic algorithm solver is run from the routine '''Gafortran'''. In this routine a set of parameter values are generated with respect to the initial parameter values, for the first run, and after that they are set to the ‘fittest’ parameters as determined by the genetic algorithm solver. The parameters supplied by Gafortran are used in the main generic glacier model, Valley, described in '''[[Model Structure]]'''.&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Optimisation_mode&amp;diff=4408</id>
		<title>Optimisation mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Optimisation_mode&amp;diff=4408"/>
		<updated>2007-11-16T17:33:13Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: /* Optimisation mode */&lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;== Optimisation mode ==&lt;br /&gt;
In this mode, the genetic algorithm solver is employed to determine the best values for the free parameters so that the difference between the modelled and the observed glacier response to climatic forcing is minimized. The cost function is determined by,&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;C = \frac{1}{n_{i}}\sum_{i} \left|{M_{i}-m_{i}}\right|+\frac{1}{n_{j}}\sum_{j} \left|{M_{j}-m_{j}}\right|+\frac{1}{n_{k}}\sum_{k} \left|{M_{k}-m_{k}}\right|+...&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
where, C is the cost function value, M is the measured value, m is the model predicted value, n is the normalisation constant, and the suffix i, j, k, refers to the different optimisation datasets. &lt;br /&gt;
&lt;br /&gt;
To run in this mode, the subroutine [[Creategadriver]] is called to create the Ga.inp file, which takes the parameters that are going to be optimised and puts them in a format that the genetic algorithm solver can read. &lt;br /&gt;
The genetic algorithm solver is run from the routine [[Gafortran]]. In this routine a set of parameter values are generated with respect to the initial parameter values, for the first run, and after that they are set to the ‘fittest’ parameters as determined by the genetic algorithm solver. The parameters supplied by Gafortran are used in the main generic glacier model, Valley, described in '''[[Model Structure]]'''.&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Optimisation_mode&amp;diff=4407</id>
		<title>Optimisation mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Optimisation_mode&amp;diff=4407"/>
		<updated>2007-11-16T17:32:58Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: /* Optimisation mode */&lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;== Optimisation mode ==&lt;br /&gt;
In this mode, the genetic algorithm solver is employed to determine the best values for the free parameters so that the difference between the modelled and the observed glacier response to climatic forcing is minimized. The cost function is determined by,&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;C = \frac{1}{n_{i}}\sum_{i} \left|{M_{i}-m_{i}}\right|+\frac{1}{n_{j}}\sum_{j} \left|{M_{j}-m_{j}}\right|+\frac{1}{n_{k}}\sum_{k} \left|{M_{k}-m_{k}}\right|+...&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
where, C is the cost function value, M is the measured value, m is the model predicted value, n is the normalisation constant, and the suffix i, j, k, refers to the different optimisation datasets. &lt;br /&gt;
&lt;br /&gt;
To run in this mode, the subroutine [[Creategadriver]] is called to create the Ga.inp file, which takes the parameters that are going to be optimised and puts them in a format that the genetic algorithm solver can read. &lt;br /&gt;
The genetic algorithm solver is run from the routine [[Gafortran]]. In this routine a set of parameter values are generated with respect to the initial parameter values, for the first run, and after that they are set to the ‘fittest’ parameters as determined by the genetic algorithm solver. The parameters supplied by Gafortran are used in the main generic glacier model, Valley, described in [[Model Structure]].&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Skadia&amp;diff=4406</id>
		<title>Skadia</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Skadia&amp;diff=4406"/>
		<updated>2007-11-16T17:32:30Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: &lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;Skadia is a 1d flow-line model that is capable of modelling glaciated systems, even with complex tributary systems. The model physics is well described in the literature (Refs.). The unique element of this model is that it uses an automated method of optimising model parameters. This methodology allows for the generic model to deal with a range of different glaciers, both in terms of the climatic and topographic setting.  &lt;br /&gt;
&lt;br /&gt;
Sections&lt;br /&gt;
#	[[Control file]]&lt;br /&gt;
#	[[Input files]]&lt;br /&gt;
#	[[Scaling]]&lt;br /&gt;
#	[[Output files]]&lt;br /&gt;
#	[[Model Structure]]&lt;br /&gt;
#	[[Optimisation]]&lt;br /&gt;
#       [[Default parameter values]]&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Skadia&amp;diff=4405</id>
		<title>Skadia</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Skadia&amp;diff=4405"/>
		<updated>2007-11-16T17:32:14Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: &lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;Skadia is a 1d flow-line model that is capable of modelling glaciated systems, even with complex tributary systems. The model physics is well described in the literature (Refs.). The unique element of this model is that it uses an automated method of optimising model parameters. This methodology allows for the generic model to deal with a range of different glaciers, both in terms of the climatic and topographic setting.  &lt;br /&gt;
&lt;br /&gt;
Sections&lt;br /&gt;
#	[[Control file]]&lt;br /&gt;
#	[[Input files]]&lt;br /&gt;
#	[[Scaling]]&lt;br /&gt;
#	[[Output files]]&lt;br /&gt;
#	[[Model structure]]&lt;br /&gt;
#	[[Optimisation]]&lt;br /&gt;
#       [[Default parameter values]]&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Optimisation_mode&amp;diff=4404</id>
		<title>Optimisation mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Optimisation_mode&amp;diff=4404"/>
		<updated>2007-11-16T17:31:48Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: /* Optimisation mode */&lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;== Optimisation mode ==&lt;br /&gt;
In this mode, the genetic algorithm solver is employed to determine the best values for the free parameters so that the difference between the modelled and the observed glacier response to climatic forcing is minimized. The cost function is determined by,&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;C = \frac{1}{n_{i}}\sum_{i} \left|{M_{i}-m_{i}}\right|+\frac{1}{n_{j}}\sum_{j} \left|{M_{j}-m_{j}}\right|+\frac{1}{n_{k}}\sum_{k} \left|{M_{k}-m_{k}}\right|+...&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
where, C is the cost function value, M is the measured value, m is the model predicted value, n is the normalisation constant, and the suffix i, j, k, refers to the different optimisation datasets. &lt;br /&gt;
&lt;br /&gt;
To run in this mode, the subroutine [[Creategadriver]] is called to create the Ga.inp file, which takes the parameters that are going to be optimised and puts them in a format that the genetic algorithm solver can read. &lt;br /&gt;
The genetic algorithm solver is run from the routine [[Gafortran]]. In this routine a set of parameter values are generated with respect to the initial parameter values, for the first run, and after that they are set to the ‘fittest’ parameters as determined by the genetic algorithm solver. The parameters supplied by Gafortran are used in the main generic glacier model, Valley, descrobed in [[model structure]].&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Optimisation_mode&amp;diff=4403</id>
		<title>Optimisation mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Optimisation_mode&amp;diff=4403"/>
		<updated>2007-11-16T17:11:40Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: /* Optimisation mode */&lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;== Optimisation mode ==&lt;br /&gt;
In this mode, the genetic algorithm solver is employed to determine the best values for the free parameters so that the difference between the modelled and the observed glacier response to climatic forcing is minimized. The cost function is determined by,&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;C = \frac{1}{n_{i}}\sum_{i} \left|{M_{i}-m_{i}}\right|+\frac{1}{n_{j}}\sum_{j} \left|{M_{j}-m_{j}}\right|+\frac{1}{n_{k}}\sum_{k} \left|{M_{k}-m_{k}}\right|+...&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
where, C is the cost function value, M is the measured value, m is the model predicted value, n is the normalisation constant, and the suffix i, j, k, refers to the different optimisation datasets. &lt;br /&gt;
&lt;br /&gt;
To run in this mode, the subroutine [[Creategadriver]] is called to create the Ga.inp file, which takes the parameters that are going to be optimised and puts them in a format that the genetic algorithm solver can read. &lt;br /&gt;
The genetic algorithm solver is run from the routine [[Gafortran]]. In this routine a set of parameter values are generated with respect to the initial parameter values, for the first run, and after that they are set to the ‘fittest’ parameters as determined by the genetic algorithm solver. The parameters supplied by Gafortran the main generic glacier model, [[Valley]], is called.&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Optimisation_mode&amp;diff=4402</id>
		<title>Optimisation mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Optimisation_mode&amp;diff=4402"/>
		<updated>2007-11-16T17:10:54Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: /* Optimisation mode */&lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;== Optimisation mode ==&lt;br /&gt;
In this mode, the genetic algorithm solver is employed to determine the best values for the free parameters so that the difference between the modelled and the observed glacier response to climatic forcing is minimized. The cost function is determined by,&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;C = \frac{1}{n_{i}}\sum_{i} \left|{M_{i}-m_{i}}\right|+\frac{1}{n_{j}}\sum_{j} \left|{M_{j}-m_{j}}\right|+\frac{1}{n_{k}}\sum_{k} \left|{M_{k}-m_{k}}\right|+...&amp;lt;/math&amp;gt;&lt;br /&gt;
&lt;br /&gt;
To run in this mode, the subroutine [[Creategadriver]] is called to create the Ga.inp file, which takes the parameters that are going to be optimised and puts them in a format that the genetic algorithm solver can read. &lt;br /&gt;
The genetic algorithm solver is run from the routine [[Gafortran]]. In this routine a set of parameter values are generated with respect to the initial parameter values, for the first run, and after that they are set to the ‘fittest’ parameters as determined by the genetic algorithm solver. The parameters supplied by Gafortran the main generic glacier model, [[Valley]], is called.&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Optimisation_mode&amp;diff=4401</id>
		<title>Optimisation mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Optimisation_mode&amp;diff=4401"/>
		<updated>2007-11-16T17:09:54Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: /* Optimisation mode */&lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;== Optimisation mode ==&lt;br /&gt;
In this mode, the genetic algorithm solver is employed to determine the best values for the free parameters so that the difference between the modelled and the observed glacier response to climatic forcing is minimized. The cost function is determined by,&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;C = \frac{1}{n_{i}}\sum_{i} \left|{M_{i}-m_{i}}\right|&amp;lt;/math&amp;gt;&lt;br /&gt;
&lt;br /&gt;
To run in this mode, the subroutine [[Creategadriver]] is called to create the Ga.inp file, which takes the parameters that are going to be optimised and puts them in a format that the genetic algorithm solver can read. &lt;br /&gt;
The genetic algorithm solver is run from the routine [[Gafortran]]. In this routine a set of parameter values are generated with respect to the initial parameter values, for the first run, and after that they are set to the ‘fittest’ parameters as determined by the genetic algorithm solver. The parameters supplied by Gafortran the main generic glacier model, [[Valley]], is called.&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Optimisation_mode&amp;diff=4400</id>
		<title>Optimisation mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Optimisation_mode&amp;diff=4400"/>
		<updated>2007-11-16T17:09:08Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: /* Optimisation mode */&lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;== Optimisation mode ==&lt;br /&gt;
In this mode, the genetic algorithm solver is employed to determine the best values for the free parameters so that the difference between the modelled and the observed glacier response to climatic forcing is minimized. The cost function is determined by,&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;C = \frac{1}{n_{i}}\sum_{i} &amp;lt;/math&amp;gt;&lt;br /&gt;
&lt;br /&gt;
To run in this mode, the subroutine [[Creategadriver]] is called to create the Ga.inp file, which takes the parameters that are going to be optimised and puts them in a format that the genetic algorithm solver can read. &lt;br /&gt;
The genetic algorithm solver is run from the routine [[Gafortran]]. In this routine a set of parameter values are generated with respect to the initial parameter values, for the first run, and after that they are set to the ‘fittest’ parameters as determined by the genetic algorithm solver. The parameters supplied by Gafortran the main generic glacier model, [[Valley]], is called.&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Optimisation_mode&amp;diff=4399</id>
		<title>Optimisation mode</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Optimisation_mode&amp;diff=4399"/>
		<updated>2007-11-16T17:08:49Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: /* Optimisation mode */&lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;== Optimisation mode ==&lt;br /&gt;
In this mode, the genetic algorithm solver is employed to determine the best values for the free parameters so that the difference between the modelled and the observed glacier response to climatic forcing is minimized. The cost function is determined by,&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;C = frac{1}{n_{i}}\sum_{i} &amp;lt;/math&amp;gt;&lt;br /&gt;
&lt;br /&gt;
To run in this mode, the subroutine [[Creategadriver]] is called to create the Ga.inp file, which takes the parameters that are going to be optimised and puts them in a format that the genetic algorithm solver can read. &lt;br /&gt;
The genetic algorithm solver is run from the routine [[Gafortran]]. In this routine a set of parameter values are generated with respect to the initial parameter values, for the first run, and after that they are set to the ‘fittest’ parameters as determined by the genetic algorithm solver. The parameters supplied by Gafortran the main generic glacier model, [[Valley]], is called.&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Model_mass_balance&amp;diff=4398</id>
		<title>Model mass balance</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Model_mass_balance&amp;diff=4398"/>
		<updated>2007-11-13T17:57:44Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: /* subroutine nexttimestep_ebm */&lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;There are two methods to calculate the mass balance within the model; &lt;br /&gt;
:*degree day method&lt;br /&gt;
:*simple energy balance model.&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
== Degree day method ==&lt;br /&gt;
The degree day method is the simplest method to calculate the ablation. In this type of model the potential ablation is determined as a linear function of the excess of daily mean air temperature above 0 °C. This potential ablation is transformed into an amount of melt through the application of degree day factors. The positive degree-day factor parameterises the physical processes that characterise the exchange of energy between the atmosphere and the surface. As a parameterisation of complex surface interchanges, degree-day factors are spatially and temporally variable. The degree day factor for snow and ice are two of the free parameters that are specified in '''[[Readoptimfile]]'''.&lt;br /&gt;
&lt;br /&gt;
=== function rainorsnw ===&lt;br /&gt;
&lt;br /&gt;
This function determines whether precipitation falls as rain or snow using air temperature(it returns amount of rainfall). If it is called with the flag not equal to 1, a simple comparison between the daily mean temperature and a prescribed rain-snow transition temperature is made: if the daily mean air temperature is above the threshold then all precipitation on that day falls as rain, otherwise it is assumed that all the precipitation is snow, so the function returns zero.&amp;lt;br&amp;gt;  &lt;br /&gt;
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 (&amp;lt;math&amp;gt;\bar{T}&amp;lt;/math&amp;gt;) with prescribed range, &amp;lt;math&amp;gt;\Delta\,T&amp;lt;/math&amp;gt;, (the range is from the mean to the peak, NOT peak to peak),&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;T(t)=\bar{T}-\Delta\,T \cos(t)&amp;lt;/math&amp;gt;&lt;br /&gt;
	 &lt;br /&gt;
which can be rearranged to give time of day &amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt; (in radians) at which the threshold &amp;lt;math&amp;gt;T^{'}&amp;lt;/math&amp;gt;  is reached, &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;t^{'}=\cos^{-1}\left({\frac{\bar{T}-T^{'}}{\Delta\,T}}\right)&amp;lt;/math&amp;gt; .&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
The proportion of the day for which temperatures exceed the threshold is therefore &amp;lt;math&amp;gt;(\pi\,-t^{'})/\pi\, &amp;lt;/math&amp;gt; .  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.&lt;br /&gt;
&lt;br /&gt;
=== function finddegdays ===&lt;br /&gt;
&lt;br /&gt;
This function finds the number of positive day degrees (''d'') 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, or none if the temperature is below zero.  &lt;br /&gt;
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 (&amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt;) with prescribed range, &amp;lt;math&amp;gt;\Delta\,T&amp;lt;/math&amp;gt;, (the range is from the mean to the peak, NOT peak to peak),&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;T(t)=\bar{T}-\Delta\,T \cos(t)&amp;lt;/math&amp;gt; &lt;br /&gt;
&lt;br /&gt;
which can be rearranged to give time of day, &amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt;, (in radians) at which the temperature equals zero, &amp;lt;br&amp;gt; &lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;t^{'}=\cos^{-1}\left({\frac{\bar{T}}{\Delta\,T}}\right)&amp;lt;/math&amp;gt; .&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
We now need to integrate ''T(t)'' from the this time to midday &amp;lt;math&amp;gt;\left({\pi\,}\right)&amp;lt;/math&amp;gt; and scale from radians into fractions of a day,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;d=\frac{1}{\pi\,} \int\limits_{t^{'}}^{\pi\,}\bar{T}-\Delta\,T\cos(t)dt&amp;lt;/math&amp;gt;&lt;br /&gt;
&lt;br /&gt;
:&amp;lt;math&amp;gt;=\frac{1}{\pi\,}\left[{\bar{T}\left({\pi\,-t^{'}}\right)+\Delta\,T\sin\left({t^{'}}\right)}\right]&amp;lt;/math&amp;gt;.&lt;br /&gt;
&lt;br /&gt;
=== function lapsecalc ===&lt;br /&gt;
&lt;br /&gt;
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 meteorological data.  &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
Although the free atmosphere temperature lapse rate is assumed to be constant at between – 6 °C to – 7 °C and this value is often used in studies investigating glacier mass balance response to climate change, the actual vertical gradient in surface temperatures is found to be spatially and temporally variable. The lapse rates are therefore opitmised in the model and so are defined in '''[[Readoptimfile]]''' . The lapse rate is applied relative to the elevation of the meteorological data, either the elevation of the meteorological station, or sea-level for ERA 40 data. All calculations are done in scaled units.&lt;br /&gt;
&lt;br /&gt;
=== function surfacemodel ===&lt;br /&gt;
&lt;br /&gt;
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.&lt;br /&gt;
# The depth of snow is increased by snowfall (total precipitation minus rainfall).&lt;br /&gt;
# The depth of super-imposed ice is increased by rainfall (assumed to refreeze).&lt;br /&gt;
# 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. &lt;br /&gt;
# 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,&lt;br /&gt;
&lt;br /&gt;
:::&amp;lt;math&amp;gt;i=w\left({s+i}\right)&amp;lt;/math&amp;gt; &lt;br /&gt;
&lt;br /&gt;
:::&amp;lt;math&amp;gt;i=\frac{ws}{1-w}&amp;lt;/math&amp;gt;&lt;br /&gt;
&amp;lt;ol&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;5&amp;quot;&amp;gt;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.&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;6&amp;quot;&amp;gt;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.&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;7&amp;quot;&amp;gt;If any melt potential remains then proceed to melt glacier ice ice (reducing potential and increasing runoff accordingly).&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;/ol&amp;gt;&lt;br /&gt;
&lt;br /&gt;
== Energy balance model ==&lt;br /&gt;
&lt;br /&gt;
=== subroutine initialstart_ebm ===&lt;br /&gt;
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.&lt;br /&gt;
&lt;br /&gt;
=== subroutine initial4day_ebm ===&lt;br /&gt;
&lt;br /&gt;
Calculates the energy-budget parameters that vary daily.  In this case, this is the solar declination (&amp;lt;math&amp;gt;\delta\,&amp;lt;/math&amp;gt;) which is then used later on with the latitude (&amp;lt;math&amp;gt;\phi\,&amp;lt;/math&amp;gt;) to calculate the solar zenith (z) and azimuth (&amp;lt;math&amp;gt;\psi\,&amp;lt;/math&amp;gt;), &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;&lt;br /&gt;
\begin{cases}&lt;br /&gt;
\delta\,= 0.006918 - 0.399912\cos\theta\,_{yr}+ 0.70257\sin\theta_{yr} -\\ 0.006758\cos2\theta_{yr} +0.000907\sin2\theta_{yr}- 0.002697\cos3\theta_{yr}+\\ 0.00148\sin3\theta_{yr} &lt;br /&gt;
\end{cases}&lt;br /&gt;
&amp;lt;/math&amp;gt;	 &lt;br /&gt;
&lt;br /&gt;
where &amp;lt;math&amp;gt;\theta_{yr}&amp;lt;/math&amp;gt; is the julian day expressed in radians. From this,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;z^{'}=&lt;br /&gt;
\begin{pmatrix}&lt;br /&gt;
  \sin\lambda\, &amp;amp; \sin\delta\, \\&lt;br /&gt;
  \cos\lambda\, &amp;amp; \cos\delta\,&lt;br /&gt;
\end{pmatrix}&lt;br /&gt;
&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
and,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;\psi^{'}=&lt;br /&gt;
\begin{pmatrix}&lt;br /&gt;
  \sin\lambda\, &amp;amp; \cos\delta\, \\&lt;br /&gt;
  \cos\lambda\, &amp;amp; \sin\delta\,&lt;br /&gt;
\end{pmatrix}&lt;br /&gt;
&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
the primes indicates that these are parameters used in the final calculation of zenith and azimuths.&lt;br /&gt;
&lt;br /&gt;
=== subroutine nexttimestep_ebm ===&lt;br /&gt;
&lt;br /&gt;
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, &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;\rho\,cd\frac{dT}{dt}=C\alpha\,^{'}I+BT_{2m}+AT &amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
::&amp;lt;math&amp;gt;I =I_{0}\left({\cos\,z\cos\phi\,+\sin\,z\sin\,\phi\,\cos\,\left({\psi\,-\beta\,}\right)}\right)&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
:&amp;lt;math&amp;gt;cos\,z=z^{'}_{1}+z^{'}_{2}\cos\,\theta_{d}&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
:&amp;lt;math&amp;gt;cos\,\psi\,=\frac{\psi\,^{'}_{1}\cos\,\theta\,_{d}-\psi\,^{'}_{2}}{\cos\,z}&amp;lt;/math&amp;gt;&lt;br /&gt;
&lt;br /&gt;
	 &lt;br /&gt;
&lt;br /&gt;
where &amp;lt;math&amp;gt;\rho\,&amp;lt;/math&amp;gt; is density, d is the depth of the thermally-active layer and c is the specific heat capacity (we currently assume snow throughout).  ''A'', ''B'' and ''C'' are free parameters. &amp;lt;math&amp;gt;\alpha\,^{'}&amp;lt;/math&amp;gt; is the co-albedo of the surface, ''I'' is the insolation at the surface; the 2-m air temperature (&amp;lt;math&amp;gt;T_{2m}&amp;lt;/math&amp;gt;) is based on the daily mean and amplitude value supplied. All temperatures are in K.&amp;lt;br&amp;gt;&lt;br /&gt;
&amp;lt;math&amp;gt;\beta\,&amp;lt;/math&amp;gt; is the surface’s aspect, &amp;lt;math&amp;gt;\phi\,&amp;lt;/math&amp;gt; is the slope (in radians)and &amp;lt;math&amp;gt;I_{0}&amp;lt;/math&amp;gt; is the insolation received at the top of the atmosphere, and is given a constant value of 1371 W&amp;lt;math&amp;gt;m^{2}&amp;lt;/math&amp;gt;. &amp;lt;math&amp;gt;\theta_{d}\,&amp;lt;/math&amp;gt; is the fraction of the day expressed in radians and relative to midday. Appropriate sets are taken to ensure that there is no insolation when the sun is below the horizon. &amp;lt;br&amp;gt;&lt;br /&gt;
 &lt;br /&gt;
We use an implicit time-marching scheme, &lt;br /&gt;
&amp;lt;math&amp;gt;T^{n}-T^{n-1}=\frac{\delta\,t}{\rho\,cd}\left({C\alpha\,^{'}I^{n}+BT_{2m}^{n}+\frac{A}{2}\left({T^{n}-T^{n-1}}\right)}\right)&amp;lt;/math&amp;gt;.&lt;br /&gt;
&lt;br /&gt;
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 in the active layer into an equivalent melt depth.  We use the trapezium rule to do this.  There are two cases.  &lt;br /&gt;
#The first is where the preceeding temperature is at the melt point also; this is easy.  &lt;br /&gt;
#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).&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Model_mass_balance&amp;diff=4397</id>
		<title>Model mass balance</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Model_mass_balance&amp;diff=4397"/>
		<updated>2007-11-13T17:56:30Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: /* subroutine nexttimestep_ebm */&lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;There are two methods to calculate the mass balance within the model; &lt;br /&gt;
:*degree day method&lt;br /&gt;
:*simple energy balance model.&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
== Degree day method ==&lt;br /&gt;
The degree day method is the simplest method to calculate the ablation. In this type of model the potential ablation is determined as a linear function of the excess of daily mean air temperature above 0 °C. This potential ablation is transformed into an amount of melt through the application of degree day factors. The positive degree-day factor parameterises the physical processes that characterise the exchange of energy between the atmosphere and the surface. As a parameterisation of complex surface interchanges, degree-day factors are spatially and temporally variable. The degree day factor for snow and ice are two of the free parameters that are specified in '''[[Readoptimfile]]'''.&lt;br /&gt;
&lt;br /&gt;
=== function rainorsnw ===&lt;br /&gt;
&lt;br /&gt;
This function determines whether precipitation falls as rain or snow using air temperature(it returns amount of rainfall). If it is called with the flag not equal to 1, a simple comparison between the daily mean temperature and a prescribed rain-snow transition temperature is made: if the daily mean air temperature is above the threshold then all precipitation on that day falls as rain, otherwise it is assumed that all the precipitation is snow, so the function returns zero.&amp;lt;br&amp;gt;  &lt;br /&gt;
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 (&amp;lt;math&amp;gt;\bar{T}&amp;lt;/math&amp;gt;) with prescribed range, &amp;lt;math&amp;gt;\Delta\,T&amp;lt;/math&amp;gt;, (the range is from the mean to the peak, NOT peak to peak),&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;T(t)=\bar{T}-\Delta\,T \cos(t)&amp;lt;/math&amp;gt;&lt;br /&gt;
	 &lt;br /&gt;
which can be rearranged to give time of day &amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt; (in radians) at which the threshold &amp;lt;math&amp;gt;T^{'}&amp;lt;/math&amp;gt;  is reached, &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;t^{'}=\cos^{-1}\left({\frac{\bar{T}-T^{'}}{\Delta\,T}}\right)&amp;lt;/math&amp;gt; .&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
The proportion of the day for which temperatures exceed the threshold is therefore &amp;lt;math&amp;gt;(\pi\,-t^{'})/\pi\, &amp;lt;/math&amp;gt; .  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.&lt;br /&gt;
&lt;br /&gt;
=== function finddegdays ===&lt;br /&gt;
&lt;br /&gt;
This function finds the number of positive day degrees (''d'') 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, or none if the temperature is below zero.  &lt;br /&gt;
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 (&amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt;) with prescribed range, &amp;lt;math&amp;gt;\Delta\,T&amp;lt;/math&amp;gt;, (the range is from the mean to the peak, NOT peak to peak),&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;T(t)=\bar{T}-\Delta\,T \cos(t)&amp;lt;/math&amp;gt; &lt;br /&gt;
&lt;br /&gt;
which can be rearranged to give time of day, &amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt;, (in radians) at which the temperature equals zero, &amp;lt;br&amp;gt; &lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;t^{'}=\cos^{-1}\left({\frac{\bar{T}}{\Delta\,T}}\right)&amp;lt;/math&amp;gt; .&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
We now need to integrate ''T(t)'' from the this time to midday &amp;lt;math&amp;gt;\left({\pi\,}\right)&amp;lt;/math&amp;gt; and scale from radians into fractions of a day,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;d=\frac{1}{\pi\,} \int\limits_{t^{'}}^{\pi\,}\bar{T}-\Delta\,T\cos(t)dt&amp;lt;/math&amp;gt;&lt;br /&gt;
&lt;br /&gt;
:&amp;lt;math&amp;gt;=\frac{1}{\pi\,}\left[{\bar{T}\left({\pi\,-t^{'}}\right)+\Delta\,T\sin\left({t^{'}}\right)}\right]&amp;lt;/math&amp;gt;.&lt;br /&gt;
&lt;br /&gt;
=== function lapsecalc ===&lt;br /&gt;
&lt;br /&gt;
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 meteorological data.  &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
Although the free atmosphere temperature lapse rate is assumed to be constant at between – 6 °C to – 7 °C and this value is often used in studies investigating glacier mass balance response to climate change, the actual vertical gradient in surface temperatures is found to be spatially and temporally variable. The lapse rates are therefore opitmised in the model and so are defined in '''[[Readoptimfile]]''' . The lapse rate is applied relative to the elevation of the meteorological data, either the elevation of the meteorological station, or sea-level for ERA 40 data. All calculations are done in scaled units.&lt;br /&gt;
&lt;br /&gt;
=== function surfacemodel ===&lt;br /&gt;
&lt;br /&gt;
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.&lt;br /&gt;
# The depth of snow is increased by snowfall (total precipitation minus rainfall).&lt;br /&gt;
# The depth of super-imposed ice is increased by rainfall (assumed to refreeze).&lt;br /&gt;
# 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. &lt;br /&gt;
# 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,&lt;br /&gt;
&lt;br /&gt;
:::&amp;lt;math&amp;gt;i=w\left({s+i}\right)&amp;lt;/math&amp;gt; &lt;br /&gt;
&lt;br /&gt;
:::&amp;lt;math&amp;gt;i=\frac{ws}{1-w}&amp;lt;/math&amp;gt;&lt;br /&gt;
&amp;lt;ol&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;5&amp;quot;&amp;gt;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.&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;6&amp;quot;&amp;gt;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.&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;7&amp;quot;&amp;gt;If any melt potential remains then proceed to melt glacier ice ice (reducing potential and increasing runoff accordingly).&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;/ol&amp;gt;&lt;br /&gt;
&lt;br /&gt;
== Energy balance model ==&lt;br /&gt;
&lt;br /&gt;
=== subroutine initialstart_ebm ===&lt;br /&gt;
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.&lt;br /&gt;
&lt;br /&gt;
=== subroutine initial4day_ebm ===&lt;br /&gt;
&lt;br /&gt;
Calculates the energy-budget parameters that vary daily.  In this case, this is the solar declination (&amp;lt;math&amp;gt;\delta\,&amp;lt;/math&amp;gt;) which is then used later on with the latitude (&amp;lt;math&amp;gt;\phi\,&amp;lt;/math&amp;gt;) to calculate the solar zenith (z) and azimuth (&amp;lt;math&amp;gt;\psi\,&amp;lt;/math&amp;gt;), &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;&lt;br /&gt;
\begin{cases}&lt;br /&gt;
\delta\,= 0.006918 - 0.399912\cos\theta\,_{yr}+ 0.70257\sin\theta_{yr} -\\ 0.006758\cos2\theta_{yr} +0.000907\sin2\theta_{yr}- 0.002697\cos3\theta_{yr}+\\ 0.00148\sin3\theta_{yr} &lt;br /&gt;
\end{cases}&lt;br /&gt;
&amp;lt;/math&amp;gt;	 &lt;br /&gt;
&lt;br /&gt;
where &amp;lt;math&amp;gt;\theta_{yr}&amp;lt;/math&amp;gt; is the julian day expressed in radians. From this,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;z^{'}=&lt;br /&gt;
\begin{pmatrix}&lt;br /&gt;
  \sin\lambda\, &amp;amp; \sin\delta\, \\&lt;br /&gt;
  \cos\lambda\, &amp;amp; \cos\delta\,&lt;br /&gt;
\end{pmatrix}&lt;br /&gt;
&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
and,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;\psi^{'}=&lt;br /&gt;
\begin{pmatrix}&lt;br /&gt;
  \sin\lambda\, &amp;amp; \cos\delta\, \\&lt;br /&gt;
  \cos\lambda\, &amp;amp; \sin\delta\,&lt;br /&gt;
\end{pmatrix}&lt;br /&gt;
&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
the primes indicates that these are parameters used in the final calculation of zenith and azimuths.&lt;br /&gt;
&lt;br /&gt;
=== subroutine nexttimestep_ebm ===&lt;br /&gt;
&lt;br /&gt;
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, &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;\rho\,cd\frac{dT}{dt}=C\alpha\,^{'}I+BT_{2m}+AT &amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
::&amp;lt;math&amp;gt;I =I_{0}\left({\cos\,z\cos\phi\,+\sin\,z\sin\,\phi\,\cos\,\left({\psi\,-\beta\,}\right)}\right)&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
:&amp;lt;math&amp;gt;cos\,z=z^{'}_{1}+z^{'}_{2}\cos\,\theta_{d}&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
:&amp;lt;math&amp;gt;cos\,\psi\,=\frac{\psi\,^{'}_{1}\cos\,\theta\,_{d}-\psi\,^{'}_{2}}{\cos\,z}&amp;lt;/math&amp;gt;&lt;br /&gt;
&lt;br /&gt;
	 &lt;br /&gt;
&lt;br /&gt;
where &amp;lt;math&amp;gt;\rho\,&amp;lt;/math&amp;gt; is density, d is the depth of the thermally-active layer and c is the specific heat capacity (we currently assume snow throughout).  ''A'', ''B'' and ''C'' are free parameters. &amp;lt;math&amp;gt;\alpha\,^{'}&amp;lt;/math&amp;gt; is the co-albedo of the surface, ''I'' is the insolation at the surface; the 2-m air temperature (&amp;lt;math&amp;gt;T_{2m}&amp;lt;/math&amp;gt;) is based on the daily mean and amplitude value supplied. All temperatures are in K.&amp;lt;br&amp;gt;&lt;br /&gt;
&amp;lt;math&amp;gt;\beta\,&amp;lt;/math&amp;gt; is the surface’s aspect, &amp;lt;math&amp;gt;\phi\,&amp;lt;/math&amp;gt; is the slope (in radians)and &amp;lt;math&amp;gt;I_{0}&amp;lt;/math&amp;gt; is the insolation received at the top of the atmosphere, and is given a constant value of 1371 W&amp;lt;math&amp;gt;m^{2}&amp;lt;/math&amp;gt;. &amp;lt;math&amp;gt;\theta_{d}\,&amp;lt;/math&amp;gt; is the fraction of the day expressed in radians and relative to midday. Appropriate sets are taken to ensure that there is no insolation when the sun is below the horizon. &amp;lt;br&amp;gt;&lt;br /&gt;
 &lt;br /&gt;
We use an implicit time-marching scheme, &lt;br /&gt;
&amp;lt;math&amp;gt;T^{n}-T^{n-1}=\frac{\delta\,t}{\rho\,cd}\left({C\alpha\,^{'}I^{n}+BT_{2m}^{n}+\frac{A}{2}\left({T^{n}-T^{n-1}}\right)}\right)&amp;lt;/math&amp;gt;&lt;br /&gt;
&lt;br /&gt;
 &lt;br /&gt;
&lt;br /&gt;
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.  &lt;br /&gt;
#The first is where the preceeding temperature is at the melt point also; this is easy.  &lt;br /&gt;
#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).&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Model_mass_balance&amp;diff=4396</id>
		<title>Model mass balance</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Model_mass_balance&amp;diff=4396"/>
		<updated>2007-11-13T17:51:10Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: /* subroutine nexttimestep_ebm */&lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;There are two methods to calculate the mass balance within the model; &lt;br /&gt;
:*degree day method&lt;br /&gt;
:*simple energy balance model.&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
== Degree day method ==&lt;br /&gt;
The degree day method is the simplest method to calculate the ablation. In this type of model the potential ablation is determined as a linear function of the excess of daily mean air temperature above 0 °C. This potential ablation is transformed into an amount of melt through the application of degree day factors. The positive degree-day factor parameterises the physical processes that characterise the exchange of energy between the atmosphere and the surface. As a parameterisation of complex surface interchanges, degree-day factors are spatially and temporally variable. The degree day factor for snow and ice are two of the free parameters that are specified in '''[[Readoptimfile]]'''.&lt;br /&gt;
&lt;br /&gt;
=== function rainorsnw ===&lt;br /&gt;
&lt;br /&gt;
This function determines whether precipitation falls as rain or snow using air temperature(it returns amount of rainfall). If it is called with the flag not equal to 1, a simple comparison between the daily mean temperature and a prescribed rain-snow transition temperature is made: if the daily mean air temperature is above the threshold then all precipitation on that day falls as rain, otherwise it is assumed that all the precipitation is snow, so the function returns zero.&amp;lt;br&amp;gt;  &lt;br /&gt;
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 (&amp;lt;math&amp;gt;\bar{T}&amp;lt;/math&amp;gt;) with prescribed range, &amp;lt;math&amp;gt;\Delta\,T&amp;lt;/math&amp;gt;, (the range is from the mean to the peak, NOT peak to peak),&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;T(t)=\bar{T}-\Delta\,T \cos(t)&amp;lt;/math&amp;gt;&lt;br /&gt;
	 &lt;br /&gt;
which can be rearranged to give time of day &amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt; (in radians) at which the threshold &amp;lt;math&amp;gt;T^{'}&amp;lt;/math&amp;gt;  is reached, &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;t^{'}=\cos^{-1}\left({\frac{\bar{T}-T^{'}}{\Delta\,T}}\right)&amp;lt;/math&amp;gt; .&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
The proportion of the day for which temperatures exceed the threshold is therefore &amp;lt;math&amp;gt;(\pi\,-t^{'})/\pi\, &amp;lt;/math&amp;gt; .  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.&lt;br /&gt;
&lt;br /&gt;
=== function finddegdays ===&lt;br /&gt;
&lt;br /&gt;
This function finds the number of positive day degrees (''d'') 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, or none if the temperature is below zero.  &lt;br /&gt;
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 (&amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt;) with prescribed range, &amp;lt;math&amp;gt;\Delta\,T&amp;lt;/math&amp;gt;, (the range is from the mean to the peak, NOT peak to peak),&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;T(t)=\bar{T}-\Delta\,T \cos(t)&amp;lt;/math&amp;gt; &lt;br /&gt;
&lt;br /&gt;
which can be rearranged to give time of day, &amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt;, (in radians) at which the temperature equals zero, &amp;lt;br&amp;gt; &lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;t^{'}=\cos^{-1}\left({\frac{\bar{T}}{\Delta\,T}}\right)&amp;lt;/math&amp;gt; .&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
We now need to integrate ''T(t)'' from the this time to midday &amp;lt;math&amp;gt;\left({\pi\,}\right)&amp;lt;/math&amp;gt; and scale from radians into fractions of a day,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;d=\frac{1}{\pi\,} \int\limits_{t^{'}}^{\pi\,}\bar{T}-\Delta\,T\cos(t)dt&amp;lt;/math&amp;gt;&lt;br /&gt;
&lt;br /&gt;
:&amp;lt;math&amp;gt;=\frac{1}{\pi\,}\left[{\bar{T}\left({\pi\,-t^{'}}\right)+\Delta\,T\sin\left({t^{'}}\right)}\right]&amp;lt;/math&amp;gt;.&lt;br /&gt;
&lt;br /&gt;
=== function lapsecalc ===&lt;br /&gt;
&lt;br /&gt;
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 meteorological data.  &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
Although the free atmosphere temperature lapse rate is assumed to be constant at between – 6 °C to – 7 °C and this value is often used in studies investigating glacier mass balance response to climate change, the actual vertical gradient in surface temperatures is found to be spatially and temporally variable. The lapse rates are therefore opitmised in the model and so are defined in '''[[Readoptimfile]]''' . The lapse rate is applied relative to the elevation of the meteorological data, either the elevation of the meteorological station, or sea-level for ERA 40 data. All calculations are done in scaled units.&lt;br /&gt;
&lt;br /&gt;
=== function surfacemodel ===&lt;br /&gt;
&lt;br /&gt;
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.&lt;br /&gt;
# The depth of snow is increased by snowfall (total precipitation minus rainfall).&lt;br /&gt;
# The depth of super-imposed ice is increased by rainfall (assumed to refreeze).&lt;br /&gt;
# 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. &lt;br /&gt;
# 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,&lt;br /&gt;
&lt;br /&gt;
:::&amp;lt;math&amp;gt;i=w\left({s+i}\right)&amp;lt;/math&amp;gt; &lt;br /&gt;
&lt;br /&gt;
:::&amp;lt;math&amp;gt;i=\frac{ws}{1-w}&amp;lt;/math&amp;gt;&lt;br /&gt;
&amp;lt;ol&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;5&amp;quot;&amp;gt;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.&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;6&amp;quot;&amp;gt;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.&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;7&amp;quot;&amp;gt;If any melt potential remains then proceed to melt glacier ice ice (reducing potential and increasing runoff accordingly).&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;/ol&amp;gt;&lt;br /&gt;
&lt;br /&gt;
== Energy balance model ==&lt;br /&gt;
&lt;br /&gt;
=== subroutine initialstart_ebm ===&lt;br /&gt;
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.&lt;br /&gt;
&lt;br /&gt;
=== subroutine initial4day_ebm ===&lt;br /&gt;
&lt;br /&gt;
Calculates the energy-budget parameters that vary daily.  In this case, this is the solar declination (&amp;lt;math&amp;gt;\delta\,&amp;lt;/math&amp;gt;) which is then used later on with the latitude (&amp;lt;math&amp;gt;\phi\,&amp;lt;/math&amp;gt;) to calculate the solar zenith (z) and azimuth (&amp;lt;math&amp;gt;\psi\,&amp;lt;/math&amp;gt;), &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;&lt;br /&gt;
\begin{cases}&lt;br /&gt;
\delta\,= 0.006918 - 0.399912\cos\theta\,_{yr}+ 0.70257\sin\theta_{yr} -\\ 0.006758\cos2\theta_{yr} +0.000907\sin2\theta_{yr}- 0.002697\cos3\theta_{yr}+\\ 0.00148\sin3\theta_{yr} &lt;br /&gt;
\end{cases}&lt;br /&gt;
&amp;lt;/math&amp;gt;	 &lt;br /&gt;
&lt;br /&gt;
where &amp;lt;math&amp;gt;\theta_{yr}&amp;lt;/math&amp;gt; is the julian day expressed in radians. From this,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;z^{'}=&lt;br /&gt;
\begin{pmatrix}&lt;br /&gt;
  \sin\lambda\, &amp;amp; \sin\delta\, \\&lt;br /&gt;
  \cos\lambda\, &amp;amp; \cos\delta\,&lt;br /&gt;
\end{pmatrix}&lt;br /&gt;
&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
and,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;\psi^{'}=&lt;br /&gt;
\begin{pmatrix}&lt;br /&gt;
  \sin\lambda\, &amp;amp; \cos\delta\, \\&lt;br /&gt;
  \cos\lambda\, &amp;amp; \sin\delta\,&lt;br /&gt;
\end{pmatrix}&lt;br /&gt;
&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
the primes indicates that these are parameters used in the final calculation of zenith and azimuths.&lt;br /&gt;
&lt;br /&gt;
=== subroutine nexttimestep_ebm ===&lt;br /&gt;
&lt;br /&gt;
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, &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;\rho\,cd\frac{dT}{dt}=C\alpha\,^{'}I+BT_{2m}+AT &amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
::&amp;lt;math&amp;gt;I =I_{0}\left({\cos\,z\cos\phi\,+\sin\,z\sin\,\phi\,\cos\,\left({\psi\,-\beta\,}\right)}\right)&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
:&amp;lt;math&amp;gt;cos\,z=z^{'}_{1}+z^{'}_{2}\cos\,\theta_{d}&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
:&amp;lt;math&amp;gt;cos\,\psi\,=\frac{\psi\,^{'}_{1}\cos\,\theta\,_{d}-\psi\,^{'}_{2}}{\cos\,z}&amp;lt;/math&amp;gt;&lt;br /&gt;
&lt;br /&gt;
	 &lt;br /&gt;
&lt;br /&gt;
where &amp;lt;math&amp;gt;\rho\,&amp;lt;/math&amp;gt; is density, d is the depth of the thermally-active layer and c is the specific heat capacity (we currently assume snow throughout).  ''A'', ''B'' and ''C'' are free parameters. &amp;lt;math&amp;gt;\alpha\,^{'}&amp;lt;/math&amp;gt; is the co-albedo of the surface, ''I'' is the insolation at the surface; the 2-m air temperature (&amp;lt;math&amp;gt;T_{2m}&amp;lt;/math&amp;gt;) is based on the daily mean and amplitude value supplied. All temperatures are in K.&amp;lt;br&amp;gt;&lt;br /&gt;
&amp;lt;math&amp;gt;\beta\,&amp;lt;/math&amp;gt; is the surface’s aspect, &amp;lt;math&amp;gt;\phi\,&amp;lt;/math&amp;gt; is the slope (in radians)and &amp;lt;math&amp;gt;I_{0}&amp;lt;/math&amp;gt; is the insolation received at the top of the atmosphere, and is given a constant value of 1371 W&amp;lt;math&amp;gt;m^{2}&amp;lt;/math&amp;gt;. .&amp;lt;math&amp;gt;\theta_{d}\,&amp;lt;/math&amp;gt; is the fraction of the day expressed in radians and relative to midday. Appropriate sets are taken to ensure that there is no insolation when the sun is below the horizon &lt;br /&gt;
We use an implicit time-marching scheme&lt;br /&gt;
&lt;br /&gt;
 &lt;br /&gt;
&lt;br /&gt;
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.  &lt;br /&gt;
#The first is where the preceeding temperature is at the melt point also; this is easy.  &lt;br /&gt;
#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).&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Model_mass_balance&amp;diff=4395</id>
		<title>Model mass balance</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Model_mass_balance&amp;diff=4395"/>
		<updated>2007-11-13T17:46:50Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: /* subroutine nexttimestep_ebm */&lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;There are two methods to calculate the mass balance within the model; &lt;br /&gt;
:*degree day method&lt;br /&gt;
:*simple energy balance model.&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
== Degree day method ==&lt;br /&gt;
The degree day method is the simplest method to calculate the ablation. In this type of model the potential ablation is determined as a linear function of the excess of daily mean air temperature above 0 °C. This potential ablation is transformed into an amount of melt through the application of degree day factors. The positive degree-day factor parameterises the physical processes that characterise the exchange of energy between the atmosphere and the surface. As a parameterisation of complex surface interchanges, degree-day factors are spatially and temporally variable. The degree day factor for snow and ice are two of the free parameters that are specified in '''[[Readoptimfile]]'''.&lt;br /&gt;
&lt;br /&gt;
=== function rainorsnw ===&lt;br /&gt;
&lt;br /&gt;
This function determines whether precipitation falls as rain or snow using air temperature(it returns amount of rainfall). If it is called with the flag not equal to 1, a simple comparison between the daily mean temperature and a prescribed rain-snow transition temperature is made: if the daily mean air temperature is above the threshold then all precipitation on that day falls as rain, otherwise it is assumed that all the precipitation is snow, so the function returns zero.&amp;lt;br&amp;gt;  &lt;br /&gt;
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 (&amp;lt;math&amp;gt;\bar{T}&amp;lt;/math&amp;gt;) with prescribed range, &amp;lt;math&amp;gt;\Delta\,T&amp;lt;/math&amp;gt;, (the range is from the mean to the peak, NOT peak to peak),&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;T(t)=\bar{T}-\Delta\,T \cos(t)&amp;lt;/math&amp;gt;&lt;br /&gt;
	 &lt;br /&gt;
which can be rearranged to give time of day &amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt; (in radians) at which the threshold &amp;lt;math&amp;gt;T^{'}&amp;lt;/math&amp;gt;  is reached, &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;t^{'}=\cos^{-1}\left({\frac{\bar{T}-T^{'}}{\Delta\,T}}\right)&amp;lt;/math&amp;gt; .&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
The proportion of the day for which temperatures exceed the threshold is therefore &amp;lt;math&amp;gt;(\pi\,-t^{'})/\pi\, &amp;lt;/math&amp;gt; .  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.&lt;br /&gt;
&lt;br /&gt;
=== function finddegdays ===&lt;br /&gt;
&lt;br /&gt;
This function finds the number of positive day degrees (''d'') 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, or none if the temperature is below zero.  &lt;br /&gt;
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 (&amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt;) with prescribed range, &amp;lt;math&amp;gt;\Delta\,T&amp;lt;/math&amp;gt;, (the range is from the mean to the peak, NOT peak to peak),&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;T(t)=\bar{T}-\Delta\,T \cos(t)&amp;lt;/math&amp;gt; &lt;br /&gt;
&lt;br /&gt;
which can be rearranged to give time of day, &amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt;, (in radians) at which the temperature equals zero, &amp;lt;br&amp;gt; &lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;t^{'}=\cos^{-1}\left({\frac{\bar{T}}{\Delta\,T}}\right)&amp;lt;/math&amp;gt; .&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
We now need to integrate ''T(t)'' from the this time to midday &amp;lt;math&amp;gt;\left({\pi\,}\right)&amp;lt;/math&amp;gt; and scale from radians into fractions of a day,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;d=\frac{1}{\pi\,} \int\limits_{t^{'}}^{\pi\,}\bar{T}-\Delta\,T\cos(t)dt&amp;lt;/math&amp;gt;&lt;br /&gt;
&lt;br /&gt;
:&amp;lt;math&amp;gt;=\frac{1}{\pi\,}\left[{\bar{T}\left({\pi\,-t^{'}}\right)+\Delta\,T\sin\left({t^{'}}\right)}\right]&amp;lt;/math&amp;gt;.&lt;br /&gt;
&lt;br /&gt;
=== function lapsecalc ===&lt;br /&gt;
&lt;br /&gt;
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 meteorological data.  &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
Although the free atmosphere temperature lapse rate is assumed to be constant at between – 6 °C to – 7 °C and this value is often used in studies investigating glacier mass balance response to climate change, the actual vertical gradient in surface temperatures is found to be spatially and temporally variable. The lapse rates are therefore opitmised in the model and so are defined in '''[[Readoptimfile]]''' . The lapse rate is applied relative to the elevation of the meteorological data, either the elevation of the meteorological station, or sea-level for ERA 40 data. All calculations are done in scaled units.&lt;br /&gt;
&lt;br /&gt;
=== function surfacemodel ===&lt;br /&gt;
&lt;br /&gt;
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.&lt;br /&gt;
# The depth of snow is increased by snowfall (total precipitation minus rainfall).&lt;br /&gt;
# The depth of super-imposed ice is increased by rainfall (assumed to refreeze).&lt;br /&gt;
# 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. &lt;br /&gt;
# 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,&lt;br /&gt;
&lt;br /&gt;
:::&amp;lt;math&amp;gt;i=w\left({s+i}\right)&amp;lt;/math&amp;gt; &lt;br /&gt;
&lt;br /&gt;
:::&amp;lt;math&amp;gt;i=\frac{ws}{1-w}&amp;lt;/math&amp;gt;&lt;br /&gt;
&amp;lt;ol&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;5&amp;quot;&amp;gt;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.&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;6&amp;quot;&amp;gt;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.&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;7&amp;quot;&amp;gt;If any melt potential remains then proceed to melt glacier ice ice (reducing potential and increasing runoff accordingly).&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;/ol&amp;gt;&lt;br /&gt;
&lt;br /&gt;
== Energy balance model ==&lt;br /&gt;
&lt;br /&gt;
=== subroutine initialstart_ebm ===&lt;br /&gt;
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.&lt;br /&gt;
&lt;br /&gt;
=== subroutine initial4day_ebm ===&lt;br /&gt;
&lt;br /&gt;
Calculates the energy-budget parameters that vary daily.  In this case, this is the solar declination (&amp;lt;math&amp;gt;\delta\,&amp;lt;/math&amp;gt;) which is then used later on with the latitude (&amp;lt;math&amp;gt;\phi\,&amp;lt;/math&amp;gt;) to calculate the solar zenith (z) and azimuth (&amp;lt;math&amp;gt;\psi\,&amp;lt;/math&amp;gt;), &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;&lt;br /&gt;
\begin{cases}&lt;br /&gt;
\delta\,= 0.006918 - 0.399912\cos\theta\,_{yr}+ 0.70257\sin\theta_{yr} -\\ 0.006758\cos2\theta_{yr} +0.000907\sin2\theta_{yr}- 0.002697\cos3\theta_{yr}+\\ 0.00148\sin3\theta_{yr} &lt;br /&gt;
\end{cases}&lt;br /&gt;
&amp;lt;/math&amp;gt;	 &lt;br /&gt;
&lt;br /&gt;
where &amp;lt;math&amp;gt;\theta_{yr}&amp;lt;/math&amp;gt; is the julian day expressed in radians. From this,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;z^{'}=&lt;br /&gt;
\begin{pmatrix}&lt;br /&gt;
  \sin\lambda\, &amp;amp; \sin\delta\, \\&lt;br /&gt;
  \cos\lambda\, &amp;amp; \cos\delta\,&lt;br /&gt;
\end{pmatrix}&lt;br /&gt;
&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
and,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;\psi^{'}=&lt;br /&gt;
\begin{pmatrix}&lt;br /&gt;
  \sin\lambda\, &amp;amp; \cos\delta\, \\&lt;br /&gt;
  \cos\lambda\, &amp;amp; \sin\delta\,&lt;br /&gt;
\end{pmatrix}&lt;br /&gt;
&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
the primes indicates that these are parameters used in the final calculation of zenith and azimuths.&lt;br /&gt;
&lt;br /&gt;
=== subroutine nexttimestep_ebm ===&lt;br /&gt;
&lt;br /&gt;
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, &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;\rho\,cd\frac{dT}{dt}=C\alpha\,^{'}I+BT_{2m}+AT &amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
::&amp;lt;math&amp;gt;I =I_{0}\left({\cos\,z\cos\phi\,+\sin\,z\sin\,\phi\,\cos\,\left({\psi\,-\beta\,}\right)}\right)&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
:&amp;lt;math&amp;gt;cos\,z=z^{'}_{1}+z^{'}_{2}\cos\,\theta_{d}&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
:&amp;lt;math&amp;gt;cos\,\psi\,=\frac{\psi\,^{'}_{1}\cos\,\theta\,_{d}-\psi\,^{'}_{2}}{\cos\,z}&amp;lt;/math&amp;gt;&lt;br /&gt;
&lt;br /&gt;
	 &lt;br /&gt;
&lt;br /&gt;
where &amp;lt;math&amp;gt;\rho\,&amp;lt;/math&amp;gt; is density, d is the depth of the thermally-active layer and c is the specific heat capacity (we currently assume snow throughout).  ''A'', ''B'' and ''C'' are free parameters. &amp;lt;math&amp;gt;\alpha\,^{'}&amp;lt;/math&amp;gt; is the co-albedo of the surface, ''I'' is the insolation at the surface; the 2-m air temperature (&amp;lt;math&amp;gt;T_{2m}&amp;lt;/math&amp;gt;) is based on the daily mean and amplitude value supplied. All temperatures are in K.&amp;lt;br&amp;gt;&lt;br /&gt;
&amp;lt;math&amp;gt;\beta\,&amp;lt;/math&amp;gt; is the surface’s aspect, &amp;lt;math&amp;gt;\phi\,&amp;lt;/math&amp;gt; is the slope (in radians)and &amp;lt;math&amp;gt;I_{0}&amp;lt;/math&amp;gt; is the insolation received at the top of the atmosphere, and is given a constant value of 1371 W&amp;lt;math&amp;gt;m^{2}&amp;lt;/math&amp;gt;.    Appropriate sets are taken to ensure that the sun does not sun when it is below the horizon.&lt;br /&gt;
and &amp;lt;math&amp;gt;\theta_{d}\,&amp;lt;/math&amp;gt; is the fraction of the day expressed in radians and relative to midday. &lt;br /&gt;
We use an implicit time-marching scheme&lt;br /&gt;
&lt;br /&gt;
 &lt;br /&gt;
&lt;br /&gt;
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.  &lt;br /&gt;
#The first is where the preceeding temperature is at the melt point also; this is easy.  &lt;br /&gt;
#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).&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Model_mass_balance&amp;diff=4394</id>
		<title>Model mass balance</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Model_mass_balance&amp;diff=4394"/>
		<updated>2007-11-13T17:46:08Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: /* subroutine nexttimestep_ebm */&lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;There are two methods to calculate the mass balance within the model; &lt;br /&gt;
:*degree day method&lt;br /&gt;
:*simple energy balance model.&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
== Degree day method ==&lt;br /&gt;
The degree day method is the simplest method to calculate the ablation. In this type of model the potential ablation is determined as a linear function of the excess of daily mean air temperature above 0 °C. This potential ablation is transformed into an amount of melt through the application of degree day factors. The positive degree-day factor parameterises the physical processes that characterise the exchange of energy between the atmosphere and the surface. As a parameterisation of complex surface interchanges, degree-day factors are spatially and temporally variable. The degree day factor for snow and ice are two of the free parameters that are specified in '''[[Readoptimfile]]'''.&lt;br /&gt;
&lt;br /&gt;
=== function rainorsnw ===&lt;br /&gt;
&lt;br /&gt;
This function determines whether precipitation falls as rain or snow using air temperature(it returns amount of rainfall). If it is called with the flag not equal to 1, a simple comparison between the daily mean temperature and a prescribed rain-snow transition temperature is made: if the daily mean air temperature is above the threshold then all precipitation on that day falls as rain, otherwise it is assumed that all the precipitation is snow, so the function returns zero.&amp;lt;br&amp;gt;  &lt;br /&gt;
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 (&amp;lt;math&amp;gt;\bar{T}&amp;lt;/math&amp;gt;) with prescribed range, &amp;lt;math&amp;gt;\Delta\,T&amp;lt;/math&amp;gt;, (the range is from the mean to the peak, NOT peak to peak),&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;T(t)=\bar{T}-\Delta\,T \cos(t)&amp;lt;/math&amp;gt;&lt;br /&gt;
	 &lt;br /&gt;
which can be rearranged to give time of day &amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt; (in radians) at which the threshold &amp;lt;math&amp;gt;T^{'}&amp;lt;/math&amp;gt;  is reached, &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;t^{'}=\cos^{-1}\left({\frac{\bar{T}-T^{'}}{\Delta\,T}}\right)&amp;lt;/math&amp;gt; .&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
The proportion of the day for which temperatures exceed the threshold is therefore &amp;lt;math&amp;gt;(\pi\,-t^{'})/\pi\, &amp;lt;/math&amp;gt; .  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.&lt;br /&gt;
&lt;br /&gt;
=== function finddegdays ===&lt;br /&gt;
&lt;br /&gt;
This function finds the number of positive day degrees (''d'') 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, or none if the temperature is below zero.  &lt;br /&gt;
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 (&amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt;) with prescribed range, &amp;lt;math&amp;gt;\Delta\,T&amp;lt;/math&amp;gt;, (the range is from the mean to the peak, NOT peak to peak),&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;T(t)=\bar{T}-\Delta\,T \cos(t)&amp;lt;/math&amp;gt; &lt;br /&gt;
&lt;br /&gt;
which can be rearranged to give time of day, &amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt;, (in radians) at which the temperature equals zero, &amp;lt;br&amp;gt; &lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;t^{'}=\cos^{-1}\left({\frac{\bar{T}}{\Delta\,T}}\right)&amp;lt;/math&amp;gt; .&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
We now need to integrate ''T(t)'' from the this time to midday &amp;lt;math&amp;gt;\left({\pi\,}\right)&amp;lt;/math&amp;gt; and scale from radians into fractions of a day,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;d=\frac{1}{\pi\,} \int\limits_{t^{'}}^{\pi\,}\bar{T}-\Delta\,T\cos(t)dt&amp;lt;/math&amp;gt;&lt;br /&gt;
&lt;br /&gt;
:&amp;lt;math&amp;gt;=\frac{1}{\pi\,}\left[{\bar{T}\left({\pi\,-t^{'}}\right)+\Delta\,T\sin\left({t^{'}}\right)}\right]&amp;lt;/math&amp;gt;.&lt;br /&gt;
&lt;br /&gt;
=== function lapsecalc ===&lt;br /&gt;
&lt;br /&gt;
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 meteorological data.  &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
Although the free atmosphere temperature lapse rate is assumed to be constant at between – 6 °C to – 7 °C and this value is often used in studies investigating glacier mass balance response to climate change, the actual vertical gradient in surface temperatures is found to be spatially and temporally variable. The lapse rates are therefore opitmised in the model and so are defined in '''[[Readoptimfile]]''' . The lapse rate is applied relative to the elevation of the meteorological data, either the elevation of the meteorological station, or sea-level for ERA 40 data. All calculations are done in scaled units.&lt;br /&gt;
&lt;br /&gt;
=== function surfacemodel ===&lt;br /&gt;
&lt;br /&gt;
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.&lt;br /&gt;
# The depth of snow is increased by snowfall (total precipitation minus rainfall).&lt;br /&gt;
# The depth of super-imposed ice is increased by rainfall (assumed to refreeze).&lt;br /&gt;
# 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. &lt;br /&gt;
# 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,&lt;br /&gt;
&lt;br /&gt;
:::&amp;lt;math&amp;gt;i=w\left({s+i}\right)&amp;lt;/math&amp;gt; &lt;br /&gt;
&lt;br /&gt;
:::&amp;lt;math&amp;gt;i=\frac{ws}{1-w}&amp;lt;/math&amp;gt;&lt;br /&gt;
&amp;lt;ol&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;5&amp;quot;&amp;gt;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.&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;6&amp;quot;&amp;gt;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.&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;7&amp;quot;&amp;gt;If any melt potential remains then proceed to melt glacier ice ice (reducing potential and increasing runoff accordingly).&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;/ol&amp;gt;&lt;br /&gt;
&lt;br /&gt;
== Energy balance model ==&lt;br /&gt;
&lt;br /&gt;
=== subroutine initialstart_ebm ===&lt;br /&gt;
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.&lt;br /&gt;
&lt;br /&gt;
=== subroutine initial4day_ebm ===&lt;br /&gt;
&lt;br /&gt;
Calculates the energy-budget parameters that vary daily.  In this case, this is the solar declination (&amp;lt;math&amp;gt;\delta\,&amp;lt;/math&amp;gt;) which is then used later on with the latitude (&amp;lt;math&amp;gt;\phi\,&amp;lt;/math&amp;gt;) to calculate the solar zenith (z) and azimuth (&amp;lt;math&amp;gt;\psi\,&amp;lt;/math&amp;gt;), &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;&lt;br /&gt;
\begin{cases}&lt;br /&gt;
\delta\,= 0.006918 - 0.399912\cos\theta\,_{yr}+ 0.70257\sin\theta_{yr} -\\ 0.006758\cos2\theta_{yr} +0.000907\sin2\theta_{yr}- 0.002697\cos3\theta_{yr}+\\ 0.00148\sin3\theta_{yr} &lt;br /&gt;
\end{cases}&lt;br /&gt;
&amp;lt;/math&amp;gt;	 &lt;br /&gt;
&lt;br /&gt;
where &amp;lt;math&amp;gt;\theta_{yr}&amp;lt;/math&amp;gt; is the julian day expressed in radians. From this,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;z^{'}=&lt;br /&gt;
\begin{pmatrix}&lt;br /&gt;
  \sin\lambda\, &amp;amp; \sin\delta\, \\&lt;br /&gt;
  \cos\lambda\, &amp;amp; \cos\delta\,&lt;br /&gt;
\end{pmatrix}&lt;br /&gt;
&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
and,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;\psi^{'}=&lt;br /&gt;
\begin{pmatrix}&lt;br /&gt;
  \sin\lambda\, &amp;amp; \cos\delta\, \\&lt;br /&gt;
  \cos\lambda\, &amp;amp; \sin\delta\,&lt;br /&gt;
\end{pmatrix}&lt;br /&gt;
&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
the primes indicates that these are parameters used in the final calculation of zenith and azimuths.&lt;br /&gt;
&lt;br /&gt;
=== subroutine nexttimestep_ebm ===&lt;br /&gt;
&lt;br /&gt;
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, &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;\rho\,cd\frac{dT}{dt}=C\alpha\,^{'}I+BT_{2m}+AT &amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
::&amp;lt;math&amp;gt;I =I_{0}\left({\cos\,z\cos\phi\,+\sin\,z\sin\,\phi\,\cos\,\left({\psi\,-\beta\,}\right)}\right)&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
:&amp;lt;math&amp;gt;cos\,z=z^{'}_{1}+z^{'}_{2}\cos\,\theta_{d}&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
:&amp;lt;math&amp;gt;cos\,\psi\,=\frac{\psi\,^{'}_{1}\cos\,\theta\,_{d}-\psi\,^{'}_{2}}{\cos\,z}&amp;lt;/math&amp;gt;&lt;br /&gt;
&lt;br /&gt;
	 &lt;br /&gt;
&lt;br /&gt;
where &amp;lt;math&amp;gt;\rho\,&amp;lt;/math&amp;gt; is density, d is the depth of the thermally-active layer and c is the specific heat capacity (we currently assume snow throughout).  ''A'', ''B'' and ''C'' are free parameters. &amp;lt;math&amp;gt;\alpha\,^{'}&amp;lt;/math&amp;gt; is the co-albedo of the surface, ''I'' is the insolation at the surface; the 2-m air temperature (&amp;lt;math&amp;gt;T_{2m}&amp;lt;/math&amp;gt;) is based on the daily mean and amplitude value supplied. All temperatures are in K.&amp;lt;br&amp;gt;&lt;br /&gt;
&amp;lt;math&amp;gt;\beta\,&amp;lt;/math&amp;gt; is the surface’s aspect, &amp;lt;math&amp;gt;\phi\,&amp;lt;/math&amp;gt; is the slope (in radians)and &amp;lt;math&amp;gt;I_{0} is the insolation received at the top of the atmosphere, and is given a constant value of 1371 W&amp;lt;math&amp;gt;m^{2}&amp;lt;/math&amp;gt;.    Appropriate sets are taken to ensure that the sun does not sun when it is below the horizon.&lt;br /&gt;
and &amp;lt;math&amp;gt;\theta_{d}\,&amp;lt;/math&amp;gt; is the fraction of the day expressed in radians and relative to midday. &lt;br /&gt;
We use an implicit time-marching scheme&lt;br /&gt;
&lt;br /&gt;
 &lt;br /&gt;
&lt;br /&gt;
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.  &lt;br /&gt;
#The first is where the preceeding temperature is at the melt point also; this is easy.  &lt;br /&gt;
#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).&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Model_mass_balance&amp;diff=4393</id>
		<title>Model mass balance</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Model_mass_balance&amp;diff=4393"/>
		<updated>2007-11-13T17:41:44Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: /* subroutine nexttimestep_ebm */&lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;There are two methods to calculate the mass balance within the model; &lt;br /&gt;
:*degree day method&lt;br /&gt;
:*simple energy balance model.&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
== Degree day method ==&lt;br /&gt;
The degree day method is the simplest method to calculate the ablation. In this type of model the potential ablation is determined as a linear function of the excess of daily mean air temperature above 0 °C. This potential ablation is transformed into an amount of melt through the application of degree day factors. The positive degree-day factor parameterises the physical processes that characterise the exchange of energy between the atmosphere and the surface. As a parameterisation of complex surface interchanges, degree-day factors are spatially and temporally variable. The degree day factor for snow and ice are two of the free parameters that are specified in '''[[Readoptimfile]]'''.&lt;br /&gt;
&lt;br /&gt;
=== function rainorsnw ===&lt;br /&gt;
&lt;br /&gt;
This function determines whether precipitation falls as rain or snow using air temperature(it returns amount of rainfall). If it is called with the flag not equal to 1, a simple comparison between the daily mean temperature and a prescribed rain-snow transition temperature is made: if the daily mean air temperature is above the threshold then all precipitation on that day falls as rain, otherwise it is assumed that all the precipitation is snow, so the function returns zero.&amp;lt;br&amp;gt;  &lt;br /&gt;
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 (&amp;lt;math&amp;gt;\bar{T}&amp;lt;/math&amp;gt;) with prescribed range, &amp;lt;math&amp;gt;\Delta\,T&amp;lt;/math&amp;gt;, (the range is from the mean to the peak, NOT peak to peak),&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;T(t)=\bar{T}-\Delta\,T \cos(t)&amp;lt;/math&amp;gt;&lt;br /&gt;
	 &lt;br /&gt;
which can be rearranged to give time of day &amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt; (in radians) at which the threshold &amp;lt;math&amp;gt;T^{'}&amp;lt;/math&amp;gt;  is reached, &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;t^{'}=\cos^{-1}\left({\frac{\bar{T}-T^{'}}{\Delta\,T}}\right)&amp;lt;/math&amp;gt; .&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
The proportion of the day for which temperatures exceed the threshold is therefore &amp;lt;math&amp;gt;(\pi\,-t^{'})/\pi\, &amp;lt;/math&amp;gt; .  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.&lt;br /&gt;
&lt;br /&gt;
=== function finddegdays ===&lt;br /&gt;
&lt;br /&gt;
This function finds the number of positive day degrees (''d'') 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, or none if the temperature is below zero.  &lt;br /&gt;
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 (&amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt;) with prescribed range, &amp;lt;math&amp;gt;\Delta\,T&amp;lt;/math&amp;gt;, (the range is from the mean to the peak, NOT peak to peak),&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;T(t)=\bar{T}-\Delta\,T \cos(t)&amp;lt;/math&amp;gt; &lt;br /&gt;
&lt;br /&gt;
which can be rearranged to give time of day, &amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt;, (in radians) at which the temperature equals zero, &amp;lt;br&amp;gt; &lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;t^{'}=\cos^{-1}\left({\frac{\bar{T}}{\Delta\,T}}\right)&amp;lt;/math&amp;gt; .&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
We now need to integrate ''T(t)'' from the this time to midday &amp;lt;math&amp;gt;\left({\pi\,}\right)&amp;lt;/math&amp;gt; and scale from radians into fractions of a day,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;d=\frac{1}{\pi\,} \int\limits_{t^{'}}^{\pi\,}\bar{T}-\Delta\,T\cos(t)dt&amp;lt;/math&amp;gt;&lt;br /&gt;
&lt;br /&gt;
:&amp;lt;math&amp;gt;=\frac{1}{\pi\,}\left[{\bar{T}\left({\pi\,-t^{'}}\right)+\Delta\,T\sin\left({t^{'}}\right)}\right]&amp;lt;/math&amp;gt;.&lt;br /&gt;
&lt;br /&gt;
=== function lapsecalc ===&lt;br /&gt;
&lt;br /&gt;
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 meteorological data.  &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
Although the free atmosphere temperature lapse rate is assumed to be constant at between – 6 °C to – 7 °C and this value is often used in studies investigating glacier mass balance response to climate change, the actual vertical gradient in surface temperatures is found to be spatially and temporally variable. The lapse rates are therefore opitmised in the model and so are defined in '''[[Readoptimfile]]''' . The lapse rate is applied relative to the elevation of the meteorological data, either the elevation of the meteorological station, or sea-level for ERA 40 data. All calculations are done in scaled units.&lt;br /&gt;
&lt;br /&gt;
=== function surfacemodel ===&lt;br /&gt;
&lt;br /&gt;
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.&lt;br /&gt;
# The depth of snow is increased by snowfall (total precipitation minus rainfall).&lt;br /&gt;
# The depth of super-imposed ice is increased by rainfall (assumed to refreeze).&lt;br /&gt;
# 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. &lt;br /&gt;
# 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,&lt;br /&gt;
&lt;br /&gt;
:::&amp;lt;math&amp;gt;i=w\left({s+i}\right)&amp;lt;/math&amp;gt; &lt;br /&gt;
&lt;br /&gt;
:::&amp;lt;math&amp;gt;i=\frac{ws}{1-w}&amp;lt;/math&amp;gt;&lt;br /&gt;
&amp;lt;ol&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;5&amp;quot;&amp;gt;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.&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;6&amp;quot;&amp;gt;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.&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;7&amp;quot;&amp;gt;If any melt potential remains then proceed to melt glacier ice ice (reducing potential and increasing runoff accordingly).&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;/ol&amp;gt;&lt;br /&gt;
&lt;br /&gt;
== Energy balance model ==&lt;br /&gt;
&lt;br /&gt;
=== subroutine initialstart_ebm ===&lt;br /&gt;
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.&lt;br /&gt;
&lt;br /&gt;
=== subroutine initial4day_ebm ===&lt;br /&gt;
&lt;br /&gt;
Calculates the energy-budget parameters that vary daily.  In this case, this is the solar declination (&amp;lt;math&amp;gt;\delta\,&amp;lt;/math&amp;gt;) which is then used later on with the latitude (&amp;lt;math&amp;gt;\phi\,&amp;lt;/math&amp;gt;) to calculate the solar zenith (z) and azimuth (&amp;lt;math&amp;gt;\psi\,&amp;lt;/math&amp;gt;), &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;&lt;br /&gt;
\begin{cases}&lt;br /&gt;
\delta\,= 0.006918 - 0.399912\cos\theta\,_{yr}+ 0.70257\sin\theta_{yr} -\\ 0.006758\cos2\theta_{yr} +0.000907\sin2\theta_{yr}- 0.002697\cos3\theta_{yr}+\\ 0.00148\sin3\theta_{yr} &lt;br /&gt;
\end{cases}&lt;br /&gt;
&amp;lt;/math&amp;gt;	 &lt;br /&gt;
&lt;br /&gt;
where &amp;lt;math&amp;gt;\theta_{yr}&amp;lt;/math&amp;gt; is the julian day expressed in radians. From this,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;z^{'}=&lt;br /&gt;
\begin{pmatrix}&lt;br /&gt;
  \sin\lambda\, &amp;amp; \sin\delta\, \\&lt;br /&gt;
  \cos\lambda\, &amp;amp; \cos\delta\,&lt;br /&gt;
\end{pmatrix}&lt;br /&gt;
&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
and,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;\psi^{'}=&lt;br /&gt;
\begin{pmatrix}&lt;br /&gt;
  \sin\lambda\, &amp;amp; \cos\delta\, \\&lt;br /&gt;
  \cos\lambda\, &amp;amp; \sin\delta\,&lt;br /&gt;
\end{pmatrix}&lt;br /&gt;
&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
the primes indicates that these are parameters used in the final calculation of zenith and azimuths.&lt;br /&gt;
&lt;br /&gt;
=== subroutine nexttimestep_ebm ===&lt;br /&gt;
&lt;br /&gt;
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, &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;\rho\,cd\frac{dT}{dt}=C\alpha\,^{'}I+BT_{2m}+AT &amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
::&amp;lt;math&amp;gt;I =I_{0}\left({\cos\,z\cos\phi\,+\sin\,z\sin\,\phi\,\cos\,\left({\psi\,-\beta\,}\right)}\right)&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
:&amp;lt;math&amp;gt;cos\,z=z^{'}_{1}+z^{'}_{2}\cos\,\theta_{d}&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
:&amp;lt;math&amp;gt;cos\,\psi\,=\frac{\psi\,^{'}_{1}\cos\,\theta\,_{d}-\psi\,^{'}_{2}}{\cos\,z}&amp;lt;/math&amp;gt;&lt;br /&gt;
&lt;br /&gt;
	 &lt;br /&gt;
&lt;br /&gt;
where &amp;lt;math&amp;gt;\rho\,&amp;lt;/math&amp;gt; is density, d is the depth of the thermally-active layer and c is the specific heat capacity (we currently assume snow throughout).  ''A'', ''B'' and ''C'' are free parameters. &amp;lt;math&amp;gt;\alpha\,^{'}&amp;lt;/math&amp;gt; is the co-albedo of the surface, ''I'' is the insolation at the surface; the 2-m air temperature (&amp;lt;math&amp;gt;T_{2m}&amp;lt;/math&amp;gt;) is based on the daily mean and amplitude value supplied. All temperatures are in K.&amp;lt;br&amp;gt;&lt;br /&gt;
xwhere &amp;lt;math&amp;gt;\beta\,&amp;lt;/math&amp;gt; is the surface’s aspect, &amp;lt;math&amp;gt;\phi\,&amp;lt;/math&amp;gt; is the slope (in radians)      Appropriate sets are taken to ensure that the sun does not sun when it is below the horizon.&lt;br /&gt;
and &amp;lt;math&amp;gt;\theta_{d}\,&amp;lt;/math&amp;gt; is the fraction of the day expressed in radians and relative to midday. &lt;br /&gt;
We use an implicit time-marching scheme&lt;br /&gt;
&lt;br /&gt;
 &lt;br /&gt;
&lt;br /&gt;
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.  &lt;br /&gt;
#The first is where the preceeding temperature is at the melt point also; this is easy.  &lt;br /&gt;
#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).&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Model_mass_balance&amp;diff=4392</id>
		<title>Model mass balance</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Model_mass_balance&amp;diff=4392"/>
		<updated>2007-11-13T17:25:26Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: /* subroutine nexttimestep_ebm */&lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;There are two methods to calculate the mass balance within the model; &lt;br /&gt;
:*degree day method&lt;br /&gt;
:*simple energy balance model.&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
== Degree day method ==&lt;br /&gt;
The degree day method is the simplest method to calculate the ablation. In this type of model the potential ablation is determined as a linear function of the excess of daily mean air temperature above 0 °C. This potential ablation is transformed into an amount of melt through the application of degree day factors. The positive degree-day factor parameterises the physical processes that characterise the exchange of energy between the atmosphere and the surface. As a parameterisation of complex surface interchanges, degree-day factors are spatially and temporally variable. The degree day factor for snow and ice are two of the free parameters that are specified in '''[[Readoptimfile]]'''.&lt;br /&gt;
&lt;br /&gt;
=== function rainorsnw ===&lt;br /&gt;
&lt;br /&gt;
This function determines whether precipitation falls as rain or snow using air temperature(it returns amount of rainfall). If it is called with the flag not equal to 1, a simple comparison between the daily mean temperature and a prescribed rain-snow transition temperature is made: if the daily mean air temperature is above the threshold then all precipitation on that day falls as rain, otherwise it is assumed that all the precipitation is snow, so the function returns zero.&amp;lt;br&amp;gt;  &lt;br /&gt;
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 (&amp;lt;math&amp;gt;\bar{T}&amp;lt;/math&amp;gt;) with prescribed range, &amp;lt;math&amp;gt;\Delta\,T&amp;lt;/math&amp;gt;, (the range is from the mean to the peak, NOT peak to peak),&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;T(t)=\bar{T}-\Delta\,T \cos(t)&amp;lt;/math&amp;gt;&lt;br /&gt;
	 &lt;br /&gt;
which can be rearranged to give time of day &amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt; (in radians) at which the threshold &amp;lt;math&amp;gt;T^{'}&amp;lt;/math&amp;gt;  is reached, &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;t^{'}=\cos^{-1}\left({\frac{\bar{T}-T^{'}}{\Delta\,T}}\right)&amp;lt;/math&amp;gt; .&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
The proportion of the day for which temperatures exceed the threshold is therefore &amp;lt;math&amp;gt;(\pi\,-t^{'})/\pi\, &amp;lt;/math&amp;gt; .  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.&lt;br /&gt;
&lt;br /&gt;
=== function finddegdays ===&lt;br /&gt;
&lt;br /&gt;
This function finds the number of positive day degrees (''d'') 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, or none if the temperature is below zero.  &lt;br /&gt;
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 (&amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt;) with prescribed range, &amp;lt;math&amp;gt;\Delta\,T&amp;lt;/math&amp;gt;, (the range is from the mean to the peak, NOT peak to peak),&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;T(t)=\bar{T}-\Delta\,T \cos(t)&amp;lt;/math&amp;gt; &lt;br /&gt;
&lt;br /&gt;
which can be rearranged to give time of day, &amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt;, (in radians) at which the temperature equals zero, &amp;lt;br&amp;gt; &lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;t^{'}=\cos^{-1}\left({\frac{\bar{T}}{\Delta\,T}}\right)&amp;lt;/math&amp;gt; .&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
We now need to integrate ''T(t)'' from the this time to midday &amp;lt;math&amp;gt;\left({\pi\,}\right)&amp;lt;/math&amp;gt; and scale from radians into fractions of a day,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;d=\frac{1}{\pi\,} \int\limits_{t^{'}}^{\pi\,}\bar{T}-\Delta\,T\cos(t)dt&amp;lt;/math&amp;gt;&lt;br /&gt;
&lt;br /&gt;
:&amp;lt;math&amp;gt;=\frac{1}{\pi\,}\left[{\bar{T}\left({\pi\,-t^{'}}\right)+\Delta\,T\sin\left({t^{'}}\right)}\right]&amp;lt;/math&amp;gt;.&lt;br /&gt;
&lt;br /&gt;
=== function lapsecalc ===&lt;br /&gt;
&lt;br /&gt;
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 meteorological data.  &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
Although the free atmosphere temperature lapse rate is assumed to be constant at between – 6 °C to – 7 °C and this value is often used in studies investigating glacier mass balance response to climate change, the actual vertical gradient in surface temperatures is found to be spatially and temporally variable. The lapse rates are therefore opitmised in the model and so are defined in '''[[Readoptimfile]]''' . The lapse rate is applied relative to the elevation of the meteorological data, either the elevation of the meteorological station, or sea-level for ERA 40 data. All calculations are done in scaled units.&lt;br /&gt;
&lt;br /&gt;
=== function surfacemodel ===&lt;br /&gt;
&lt;br /&gt;
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.&lt;br /&gt;
# The depth of snow is increased by snowfall (total precipitation minus rainfall).&lt;br /&gt;
# The depth of super-imposed ice is increased by rainfall (assumed to refreeze).&lt;br /&gt;
# 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. &lt;br /&gt;
# 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,&lt;br /&gt;
&lt;br /&gt;
:::&amp;lt;math&amp;gt;i=w\left({s+i}\right)&amp;lt;/math&amp;gt; &lt;br /&gt;
&lt;br /&gt;
:::&amp;lt;math&amp;gt;i=\frac{ws}{1-w}&amp;lt;/math&amp;gt;&lt;br /&gt;
&amp;lt;ol&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;5&amp;quot;&amp;gt;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.&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;6&amp;quot;&amp;gt;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.&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;7&amp;quot;&amp;gt;If any melt potential remains then proceed to melt glacier ice ice (reducing potential and increasing runoff accordingly).&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;/ol&amp;gt;&lt;br /&gt;
&lt;br /&gt;
== Energy balance model ==&lt;br /&gt;
&lt;br /&gt;
=== subroutine initialstart_ebm ===&lt;br /&gt;
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.&lt;br /&gt;
&lt;br /&gt;
=== subroutine initial4day_ebm ===&lt;br /&gt;
&lt;br /&gt;
Calculates the energy-budget parameters that vary daily.  In this case, this is the solar declination (&amp;lt;math&amp;gt;\delta\,&amp;lt;/math&amp;gt;) which is then used later on with the latitude (&amp;lt;math&amp;gt;\phi\,&amp;lt;/math&amp;gt;) to calculate the solar zenith (z) and azimuth (&amp;lt;math&amp;gt;\psi\,&amp;lt;/math&amp;gt;), &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;&lt;br /&gt;
\begin{cases}&lt;br /&gt;
\delta\,= 0.006918 - 0.399912\cos\theta\,_{yr}+ 0.70257\sin\theta_{yr} -\\ 0.006758\cos2\theta_{yr} +0.000907\sin2\theta_{yr}- 0.002697\cos3\theta_{yr}+\\ 0.00148\sin3\theta_{yr} &lt;br /&gt;
\end{cases}&lt;br /&gt;
&amp;lt;/math&amp;gt;	 &lt;br /&gt;
&lt;br /&gt;
where &amp;lt;math&amp;gt;\theta_{yr}&amp;lt;/math&amp;gt; is the julian day expressed in radians. From this,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;z^{'}=&lt;br /&gt;
\begin{pmatrix}&lt;br /&gt;
  \sin\lambda\, &amp;amp; \sin\delta\, \\&lt;br /&gt;
  \cos\lambda\, &amp;amp; \cos\delta\,&lt;br /&gt;
\end{pmatrix}&lt;br /&gt;
&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
and,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;\psi^{'}=&lt;br /&gt;
\begin{pmatrix}&lt;br /&gt;
  \sin\lambda\, &amp;amp; \cos\delta\, \\&lt;br /&gt;
  \cos\lambda\, &amp;amp; \sin\delta\,&lt;br /&gt;
\end{pmatrix}&lt;br /&gt;
&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
the primes indicates that these are parameters used in the final calculation of zenith and azimuths.&lt;br /&gt;
&lt;br /&gt;
=== subroutine nexttimestep_ebm ===&lt;br /&gt;
&lt;br /&gt;
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, &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;\rho\,cd\frac{dT}{dt}=C\alpha\,^{'}I+BT_{2m}+AT &amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
::&amp;lt;math&amp;gt;I =I_{0}\left({\cos\,z\cos\phi\,+\sin\,z\sin\,\phi\,\cos\,\left({\psi\,-\beta\,}\right)}\right)&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
:&amp;lt;math&amp;gt;cos\,z=z^{'}_{1}+z^{'}_{2}\cos\,\theta_{d}&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
:&amp;lt;math&amp;gt;cos\,\psi\,=\frac{\psi\,^{'}_{1}\cos\,\theta\,_{d}-\psi\,^{'}_{2}}{\cos\,z}&amp;lt;/math&amp;gt;&lt;br /&gt;
&lt;br /&gt;
	 &lt;br /&gt;
&lt;br /&gt;
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.&lt;br /&gt;
&lt;br /&gt;
We use an implicit time-marching scheme&lt;br /&gt;
&lt;br /&gt;
 &lt;br /&gt;
&lt;br /&gt;
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.  &lt;br /&gt;
#The first is where the preceeding temperature is at the melt point also; this is easy.  &lt;br /&gt;
#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).&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Model_mass_balance&amp;diff=4391</id>
		<title>Model mass balance</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Model_mass_balance&amp;diff=4391"/>
		<updated>2007-11-13T17:25:08Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: /* subroutine nexttimestep_ebm */&lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;There are two methods to calculate the mass balance within the model; &lt;br /&gt;
:*degree day method&lt;br /&gt;
:*simple energy balance model.&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
== Degree day method ==&lt;br /&gt;
The degree day method is the simplest method to calculate the ablation. In this type of model the potential ablation is determined as a linear function of the excess of daily mean air temperature above 0 °C. This potential ablation is transformed into an amount of melt through the application of degree day factors. The positive degree-day factor parameterises the physical processes that characterise the exchange of energy between the atmosphere and the surface. As a parameterisation of complex surface interchanges, degree-day factors are spatially and temporally variable. The degree day factor for snow and ice are two of the free parameters that are specified in '''[[Readoptimfile]]'''.&lt;br /&gt;
&lt;br /&gt;
=== function rainorsnw ===&lt;br /&gt;
&lt;br /&gt;
This function determines whether precipitation falls as rain or snow using air temperature(it returns amount of rainfall). If it is called with the flag not equal to 1, a simple comparison between the daily mean temperature and a prescribed rain-snow transition temperature is made: if the daily mean air temperature is above the threshold then all precipitation on that day falls as rain, otherwise it is assumed that all the precipitation is snow, so the function returns zero.&amp;lt;br&amp;gt;  &lt;br /&gt;
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 (&amp;lt;math&amp;gt;\bar{T}&amp;lt;/math&amp;gt;) with prescribed range, &amp;lt;math&amp;gt;\Delta\,T&amp;lt;/math&amp;gt;, (the range is from the mean to the peak, NOT peak to peak),&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;T(t)=\bar{T}-\Delta\,T \cos(t)&amp;lt;/math&amp;gt;&lt;br /&gt;
	 &lt;br /&gt;
which can be rearranged to give time of day &amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt; (in radians) at which the threshold &amp;lt;math&amp;gt;T^{'}&amp;lt;/math&amp;gt;  is reached, &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;t^{'}=\cos^{-1}\left({\frac{\bar{T}-T^{'}}{\Delta\,T}}\right)&amp;lt;/math&amp;gt; .&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
The proportion of the day for which temperatures exceed the threshold is therefore &amp;lt;math&amp;gt;(\pi\,-t^{'})/\pi\, &amp;lt;/math&amp;gt; .  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.&lt;br /&gt;
&lt;br /&gt;
=== function finddegdays ===&lt;br /&gt;
&lt;br /&gt;
This function finds the number of positive day degrees (''d'') 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, or none if the temperature is below zero.  &lt;br /&gt;
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 (&amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt;) with prescribed range, &amp;lt;math&amp;gt;\Delta\,T&amp;lt;/math&amp;gt;, (the range is from the mean to the peak, NOT peak to peak),&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;T(t)=\bar{T}-\Delta\,T \cos(t)&amp;lt;/math&amp;gt; &lt;br /&gt;
&lt;br /&gt;
which can be rearranged to give time of day, &amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt;, (in radians) at which the temperature equals zero, &amp;lt;br&amp;gt; &lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;t^{'}=\cos^{-1}\left({\frac{\bar{T}}{\Delta\,T}}\right)&amp;lt;/math&amp;gt; .&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
We now need to integrate ''T(t)'' from the this time to midday &amp;lt;math&amp;gt;\left({\pi\,}\right)&amp;lt;/math&amp;gt; and scale from radians into fractions of a day,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;d=\frac{1}{\pi\,} \int\limits_{t^{'}}^{\pi\,}\bar{T}-\Delta\,T\cos(t)dt&amp;lt;/math&amp;gt;&lt;br /&gt;
&lt;br /&gt;
:&amp;lt;math&amp;gt;=\frac{1}{\pi\,}\left[{\bar{T}\left({\pi\,-t^{'}}\right)+\Delta\,T\sin\left({t^{'}}\right)}\right]&amp;lt;/math&amp;gt;.&lt;br /&gt;
&lt;br /&gt;
=== function lapsecalc ===&lt;br /&gt;
&lt;br /&gt;
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 meteorological data.  &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
Although the free atmosphere temperature lapse rate is assumed to be constant at between – 6 °C to – 7 °C and this value is often used in studies investigating glacier mass balance response to climate change, the actual vertical gradient in surface temperatures is found to be spatially and temporally variable. The lapse rates are therefore opitmised in the model and so are defined in '''[[Readoptimfile]]''' . The lapse rate is applied relative to the elevation of the meteorological data, either the elevation of the meteorological station, or sea-level for ERA 40 data. All calculations are done in scaled units.&lt;br /&gt;
&lt;br /&gt;
=== function surfacemodel ===&lt;br /&gt;
&lt;br /&gt;
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.&lt;br /&gt;
# The depth of snow is increased by snowfall (total precipitation minus rainfall).&lt;br /&gt;
# The depth of super-imposed ice is increased by rainfall (assumed to refreeze).&lt;br /&gt;
# 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. &lt;br /&gt;
# 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,&lt;br /&gt;
&lt;br /&gt;
:::&amp;lt;math&amp;gt;i=w\left({s+i}\right)&amp;lt;/math&amp;gt; &lt;br /&gt;
&lt;br /&gt;
:::&amp;lt;math&amp;gt;i=\frac{ws}{1-w}&amp;lt;/math&amp;gt;&lt;br /&gt;
&amp;lt;ol&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;5&amp;quot;&amp;gt;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.&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;6&amp;quot;&amp;gt;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.&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;7&amp;quot;&amp;gt;If any melt potential remains then proceed to melt glacier ice ice (reducing potential and increasing runoff accordingly).&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;/ol&amp;gt;&lt;br /&gt;
&lt;br /&gt;
== Energy balance model ==&lt;br /&gt;
&lt;br /&gt;
=== subroutine initialstart_ebm ===&lt;br /&gt;
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.&lt;br /&gt;
&lt;br /&gt;
=== subroutine initial4day_ebm ===&lt;br /&gt;
&lt;br /&gt;
Calculates the energy-budget parameters that vary daily.  In this case, this is the solar declination (&amp;lt;math&amp;gt;\delta\,&amp;lt;/math&amp;gt;) which is then used later on with the latitude (&amp;lt;math&amp;gt;\phi\,&amp;lt;/math&amp;gt;) to calculate the solar zenith (z) and azimuth (&amp;lt;math&amp;gt;\psi\,&amp;lt;/math&amp;gt;), &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;&lt;br /&gt;
\begin{cases}&lt;br /&gt;
\delta\,= 0.006918 - 0.399912\cos\theta\,_{yr}+ 0.70257\sin\theta_{yr} -\\ 0.006758\cos2\theta_{yr} +0.000907\sin2\theta_{yr}- 0.002697\cos3\theta_{yr}+\\ 0.00148\sin3\theta_{yr} &lt;br /&gt;
\end{cases}&lt;br /&gt;
&amp;lt;/math&amp;gt;	 &lt;br /&gt;
&lt;br /&gt;
where &amp;lt;math&amp;gt;\theta_{yr}&amp;lt;/math&amp;gt; is the julian day expressed in radians. From this,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;z^{'}=&lt;br /&gt;
\begin{pmatrix}&lt;br /&gt;
  \sin\lambda\, &amp;amp; \sin\delta\, \\&lt;br /&gt;
  \cos\lambda\, &amp;amp; \cos\delta\,&lt;br /&gt;
\end{pmatrix}&lt;br /&gt;
&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
and,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;\psi^{'}=&lt;br /&gt;
\begin{pmatrix}&lt;br /&gt;
  \sin\lambda\, &amp;amp; \cos\delta\, \\&lt;br /&gt;
  \cos\lambda\, &amp;amp; \sin\delta\,&lt;br /&gt;
\end{pmatrix}&lt;br /&gt;
&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
the primes indicates that these are parameters used in the final calculation of zenith and azimuths.&lt;br /&gt;
&lt;br /&gt;
=== subroutine nexttimestep_ebm ===&lt;br /&gt;
&lt;br /&gt;
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, &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;\rho\,cd\frac{dT}{dt}=C\alpha\,^{'}I+BT_{2m}+AT &amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
::&amp;lt;math&amp;gt;I =I_{0}\left({\cos\,z\cos\phi\,+\sin\,z\sin\,\phi\,\cos\,\left({\psi\,-\beta\,}\right)}\right)&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
:&amp;lt;math&amp;gt;cos\,z=z^{'}_{1}+z^{'}_{2}\cos\,\theta_{d}&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
:&amp;lt;math&amp;gt;cos\,\psi\,=\frac{\psi\,^{'}_{1}\cos\,\theta\,_{d}-\psi\,^{'}_{2}}{\cos\,z}&amp;lt;/math&amp;gt;&lt;br /&gt;
&lt;br /&gt;
	 &lt;br /&gt;
&lt;br /&gt;
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.&lt;br /&gt;
&lt;br /&gt;
We use an implicit time-marching scheme&lt;br /&gt;
&lt;br /&gt;
 &lt;br /&gt;
&lt;br /&gt;
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.  &lt;br /&gt;
#The first is where the preceeding temperature is at the melt point also; this is easy.  &lt;br /&gt;
#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).&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
	<entry>
		<id>https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Model_mass_balance&amp;diff=4390</id>
		<title>Model mass balance</title>
		<link rel="alternate" type="text/html" href="https://source.geography.bristol.ac.uk/mediawiki/index.php?title=Model_mass_balance&amp;diff=4390"/>
		<updated>2007-11-13T17:21:50Z</updated>

		<summary type="html">&lt;p&gt;RosDeath: /* subroutine initial4day_ebm */&lt;/p&gt;
&lt;hr /&gt;
&lt;div&gt;There are two methods to calculate the mass balance within the model; &lt;br /&gt;
:*degree day method&lt;br /&gt;
:*simple energy balance model.&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
== Degree day method ==&lt;br /&gt;
The degree day method is the simplest method to calculate the ablation. In this type of model the potential ablation is determined as a linear function of the excess of daily mean air temperature above 0 °C. This potential ablation is transformed into an amount of melt through the application of degree day factors. The positive degree-day factor parameterises the physical processes that characterise the exchange of energy between the atmosphere and the surface. As a parameterisation of complex surface interchanges, degree-day factors are spatially and temporally variable. The degree day factor for snow and ice are two of the free parameters that are specified in '''[[Readoptimfile]]'''.&lt;br /&gt;
&lt;br /&gt;
=== function rainorsnw ===&lt;br /&gt;
&lt;br /&gt;
This function determines whether precipitation falls as rain or snow using air temperature(it returns amount of rainfall). If it is called with the flag not equal to 1, a simple comparison between the daily mean temperature and a prescribed rain-snow transition temperature is made: if the daily mean air temperature is above the threshold then all precipitation on that day falls as rain, otherwise it is assumed that all the precipitation is snow, so the function returns zero.&amp;lt;br&amp;gt;  &lt;br /&gt;
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 (&amp;lt;math&amp;gt;\bar{T}&amp;lt;/math&amp;gt;) with prescribed range, &amp;lt;math&amp;gt;\Delta\,T&amp;lt;/math&amp;gt;, (the range is from the mean to the peak, NOT peak to peak),&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;T(t)=\bar{T}-\Delta\,T \cos(t)&amp;lt;/math&amp;gt;&lt;br /&gt;
	 &lt;br /&gt;
which can be rearranged to give time of day &amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt; (in radians) at which the threshold &amp;lt;math&amp;gt;T^{'}&amp;lt;/math&amp;gt;  is reached, &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;t^{'}=\cos^{-1}\left({\frac{\bar{T}-T^{'}}{\Delta\,T}}\right)&amp;lt;/math&amp;gt; .&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
The proportion of the day for which temperatures exceed the threshold is therefore &amp;lt;math&amp;gt;(\pi\,-t^{'})/\pi\, &amp;lt;/math&amp;gt; .  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.&lt;br /&gt;
&lt;br /&gt;
=== function finddegdays ===&lt;br /&gt;
&lt;br /&gt;
This function finds the number of positive day degrees (''d'') 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, or none if the temperature is below zero.  &lt;br /&gt;
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 (&amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt;) with prescribed range, &amp;lt;math&amp;gt;\Delta\,T&amp;lt;/math&amp;gt;, (the range is from the mean to the peak, NOT peak to peak),&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;T(t)=\bar{T}-\Delta\,T \cos(t)&amp;lt;/math&amp;gt; &lt;br /&gt;
&lt;br /&gt;
which can be rearranged to give time of day, &amp;lt;math&amp;gt;t^{'}&amp;lt;/math&amp;gt;, (in radians) at which the temperature equals zero, &amp;lt;br&amp;gt; &lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;t^{'}=\cos^{-1}\left({\frac{\bar{T}}{\Delta\,T}}\right)&amp;lt;/math&amp;gt; .&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
We now need to integrate ''T(t)'' from the this time to midday &amp;lt;math&amp;gt;\left({\pi\,}\right)&amp;lt;/math&amp;gt; and scale from radians into fractions of a day,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;d=\frac{1}{\pi\,} \int\limits_{t^{'}}^{\pi\,}\bar{T}-\Delta\,T\cos(t)dt&amp;lt;/math&amp;gt;&lt;br /&gt;
&lt;br /&gt;
:&amp;lt;math&amp;gt;=\frac{1}{\pi\,}\left[{\bar{T}\left({\pi\,-t^{'}}\right)+\Delta\,T\sin\left({t^{'}}\right)}\right]&amp;lt;/math&amp;gt;.&lt;br /&gt;
&lt;br /&gt;
=== function lapsecalc ===&lt;br /&gt;
&lt;br /&gt;
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 meteorological data.  &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
Although the free atmosphere temperature lapse rate is assumed to be constant at between – 6 °C to – 7 °C and this value is often used in studies investigating glacier mass balance response to climate change, the actual vertical gradient in surface temperatures is found to be spatially and temporally variable. The lapse rates are therefore opitmised in the model and so are defined in '''[[Readoptimfile]]''' . The lapse rate is applied relative to the elevation of the meteorological data, either the elevation of the meteorological station, or sea-level for ERA 40 data. All calculations are done in scaled units.&lt;br /&gt;
&lt;br /&gt;
=== function surfacemodel ===&lt;br /&gt;
&lt;br /&gt;
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.&lt;br /&gt;
# The depth of snow is increased by snowfall (total precipitation minus rainfall).&lt;br /&gt;
# The depth of super-imposed ice is increased by rainfall (assumed to refreeze).&lt;br /&gt;
# 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. &lt;br /&gt;
# 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,&lt;br /&gt;
&lt;br /&gt;
:::&amp;lt;math&amp;gt;i=w\left({s+i}\right)&amp;lt;/math&amp;gt; &lt;br /&gt;
&lt;br /&gt;
:::&amp;lt;math&amp;gt;i=\frac{ws}{1-w}&amp;lt;/math&amp;gt;&lt;br /&gt;
&amp;lt;ol&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;5&amp;quot;&amp;gt;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.&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;6&amp;quot;&amp;gt;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.&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;li value=&amp;quot;7&amp;quot;&amp;gt;If any melt potential remains then proceed to melt glacier ice ice (reducing potential and increasing runoff accordingly).&amp;lt;/li&amp;gt;&lt;br /&gt;
&amp;lt;/ol&amp;gt;&lt;br /&gt;
&lt;br /&gt;
== Energy balance model ==&lt;br /&gt;
&lt;br /&gt;
=== subroutine initialstart_ebm ===&lt;br /&gt;
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.&lt;br /&gt;
&lt;br /&gt;
=== subroutine initial4day_ebm ===&lt;br /&gt;
&lt;br /&gt;
Calculates the energy-budget parameters that vary daily.  In this case, this is the solar declination (&amp;lt;math&amp;gt;\delta\,&amp;lt;/math&amp;gt;) which is then used later on with the latitude (&amp;lt;math&amp;gt;\phi\,&amp;lt;/math&amp;gt;) to calculate the solar zenith (z) and azimuth (&amp;lt;math&amp;gt;\psi\,&amp;lt;/math&amp;gt;), &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;&lt;br /&gt;
\begin{cases}&lt;br /&gt;
\delta\,= 0.006918 - 0.399912\cos\theta\,_{yr}+ 0.70257\sin\theta_{yr} -\\ 0.006758\cos2\theta_{yr} +0.000907\sin2\theta_{yr}- 0.002697\cos3\theta_{yr}+\\ 0.00148\sin3\theta_{yr} &lt;br /&gt;
\end{cases}&lt;br /&gt;
&amp;lt;/math&amp;gt;	 &lt;br /&gt;
&lt;br /&gt;
where &amp;lt;math&amp;gt;\theta_{yr}&amp;lt;/math&amp;gt; is the julian day expressed in radians. From this,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;z^{'}=&lt;br /&gt;
\begin{pmatrix}&lt;br /&gt;
  \sin\lambda\, &amp;amp; \sin\delta\, \\&lt;br /&gt;
  \cos\lambda\, &amp;amp; \cos\delta\,&lt;br /&gt;
\end{pmatrix}&lt;br /&gt;
&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
and,&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;\psi^{'}=&lt;br /&gt;
\begin{pmatrix}&lt;br /&gt;
  \sin\lambda\, &amp;amp; \cos\delta\, \\&lt;br /&gt;
  \cos\lambda\, &amp;amp; \sin\delta\,&lt;br /&gt;
\end{pmatrix}&lt;br /&gt;
&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
the primes indicates that these are parameters used in the final calculation of zenith and azimuths.&lt;br /&gt;
&lt;br /&gt;
=== subroutine nexttimestep_ebm ===&lt;br /&gt;
&lt;br /&gt;
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, &amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;math&amp;gt;\rho\,cd\frac{dT}{dt}=C\alpha\,^{'}I+BT_{2m}+AT &amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
::&amp;lt;math&amp;gt;I =I_{0}\left({\cos\,z\cos\phi\,+\sin\,z\sin\,\phi\,\cos\,\left({\psi\,-\beta\,}\right)}\right)&amp;lt;/math&amp;gt;&amp;lt;br&amp;gt;&lt;br /&gt;
&lt;br /&gt;
:&amp;lt;math&amp;gt;cos\,z=z^{'}_{1}+z^{'}_{2}\cos\,\theta_{d}&amp;lt;/math&amp;gt;&lt;br /&gt;
&lt;br /&gt;
	 &lt;br /&gt;
&lt;br /&gt;
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.&lt;br /&gt;
&lt;br /&gt;
We use an implicit time-marching scheme&lt;br /&gt;
&lt;br /&gt;
 &lt;br /&gt;
&lt;br /&gt;
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.  &lt;br /&gt;
#The first is where the preceeding temperature is at the melt point also; this is easy.  &lt;br /&gt;
#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).&lt;/div&gt;</summary>
		<author><name>RosDeath</name></author>
	</entry>
</feed>