implementation of Robinson, Walter A. (2001) Modeling Dynamic Climate Systems. (section 2.4)
<?xml version="1.0" encoding="UTF-8"?>
<!--Xholon Workbook http://www.primordion.com/Xholon/wb/ (C) Ken Webb Tue May 08 2012 08:29:51 GMT-0400 (EDT)-->
<XholonWorkbook>
<Notes><![CDATA[
Xholon
------
Title: Radiative-Convective Model of the Atmosphere
Description: implementation of Robinson, Walter A. (2001) Modeling Dynamic Climate Systems. (section 2.4)
Url: http://www.primordion.com/Xholon/wb/
InternalName:
YoutubeId:
Keywords:
My Notes
--------
How to run the app contained in this workbook ::
Click the Run button above.
Click the Step button in the overlay, once.
The overlay will display:
the n layers of the atmosphere (5 by default),
plus a layer for the top of the atmosphere (tpf),
and another layer for the surface of the Earth (srf),
and symbols showing the down and up flows of energy between layers.
Move the mouse over the graphic to see current information about each layer.
How to change the number of layers in the Atmosphere ::
Scroll down to the TheSystem editor.
Change the Layer multiplicity to some other value, such as 10 or 100.
Click the Run button (or Refresh in the overlay), and then Step.
|This Xholon Workbook is not yet complete|
Robinson discusses the flows of visible and infrared radiation through the atmosphere. He describes it as "hugely complicated". The scientific field is called radiative transfer.
He uses STELLA as a modeling tool. Figure 2.9 in the book looks much more complicated than I think it should be, so I'm going to re-implement it using Xholon, in this workbook. The Robinson figure and model look complicated because in STELLA each layer of the Atmosphere must be separately drawn, and each duplicate part separately named. With Xholon I can just specify a multiplicity. Robinson uses 5 vertical layers. Once I define the characteristics of a AtmosphereLayer and its behavior, I should be able to just specify ::
<AtmosphereLayer multiplicity="5"/>
or even ::
<AtmosphereLayer multiplicity="100"/>
and it should work.
In the model, there are three flows of radiative energy, and one convective flow ::
incoming solar radiation
down-flowing infrared radiation
up-flowing infrared radiation
convection
A convection oven also has these two basic types of flows - radiant energy from hot elements at the top and bottom, and convection from a small fan.
Convection
----------
For convection, Robinson introduces a z parameter for each layer, which he calls a "convective adjustment". z for layers 1 to 5 in STELLA notation are ::
z_1 = LOGN(1000 / 900) * (T_1 + T_0) / 2 * 29.3 / 1000
z_2 = z_1 + LOGN(900 / 700) * (T_2 + T_1) / 2 * 29.3 / 1000
z_3 = z_2 + LOGN(700 / 500) * (T_3 + T_2) / 2 * 29.3 / 1000
z_4 = z_3 + LOGN(500 / 300) * (T_4 + T_3) / 2 * 29.3 / 1000
z_5 = z_4 + LOGN(300 / 100) * (T_5 + T_4) / 2 * 29.3 / 1000
Robinson doesn't explain where these numbers come from, except to say that they produce results "much closer to what is observed".
The JavaScript implementation of these calculations is included in the "JavaScript code" editor below. To try this out, click the Run button above, select the lines of code from the editor (not including the initial and final "script" lines), copy-and-paste them (or drag them) into the textarea in the runtime overlay, and click Submit. The results are ::
z_1 = 0.8821128480167159
z_2 = 2.900819868053781
z_3 = 5.485261435180134
z_4 = 9.248761555891715
z_5 = 16.919481291707307
These are the values (rounded off) that used later in this workbook in the Atmospherebehavior editor.
Results
-------
How to check the results ::
Open the "Workbook Scripts" workbook
http://www.primordion.com/Xholon/wb/openwb.php?q=2567985&f=gist.github.com/raw/
Copy the contents of the "Show the contents of the entire CSH tree (webEdition)" script
Paste that into the textarea in the runtime overlay on this page
Click the Submit button
How to display these results organized as a mind map ::
Copy the results from the runtime textarea, starting with the "TheSystem_0" line
Go to the Mappio site at http://mappio.com/
Click the Upload MindMap button
Log in as a guest, or using your own account if you have one
Paste the copied results into the Comments area
Click the Save button
Mappio creates and displays a JPEG image that you can save
A sample saved mind map, of an early and not-quite-correct version of the model, can be viewed at either of these URLs
http://mappio.com/mindmap/guest/new-mindmap-5-6-2012-12-44-28-pm
http://www.primordion.com/Xholon/wb/images/new-mindmap-5-6-2012-12-44-28-pm-Large.jpg
Graphics
--------
- use Raphael ?
- draw a stack of rectangles, with Tpf at top, n layers in the middle, and Srf at the bottom
- draw directed triangles between layers to show radiation flux
- tooltips that show the latest values
Links
-----
Wikipedia has some potentially useful links, but to me they mostly re-emphasize Robinson's statement that radiative transfer is "hugely complicated" ::
http://en.wikipedia.org/wiki/Radiative_transfer
http://en.wikipedia.org/wiki/Atmospheric_radiative_transfer_codes
http://en.wikipedia.org/wiki/Parametrization_(climate)
http://en.wikipedia.org/wiki/Radiative_equilibrium
http://en.wikipedia.org/wiki/Kirchhoff's_law_of_thermal_radiation
]]></Notes>
<script implName="lang:python:inline:"><![CDATA[
print "Radiation"
]]></script>
<script implName="lang:javascript:inline:"><![CDATA[
print("\n Calculation of z parameters for convection:");
// initial temperatures T in degrees Kelvin (K)
var T_0 = 290.19; // temperature at the base of layer 1; this is the temperature at the Earth's surface
var T_1 = 281.3; // temperature at the top of layer 1
var T_2 = 267.0;
var T_3 = 257.3;
var T_4 = 245.6;
var T_5 = 231.0; // temperature at the top of layer 5; this is the temperature at the top of Earth's atmosphere
// calculate z for each layer
var z_1 = Math.log(1000 / 900) * (T_1 + T_0) / 2 * 29.3 / 1000;
var z_2 = z_1 + Math.log(900 / 700) * (T_2 + T_1) / 2 * 29.3 / 1000;
var z_3 = z_2 + Math.log(700 / 500) * (T_3 + T_2) / 2 * 29.3 / 1000;
var z_4 = z_3 + Math.log(500 / 300) * (T_4 + T_3) / 2 * 29.3 / 1000;
var z_5 = z_4 + Math.log(300 / 100) * (T_5 + T_4) / 2 * 29.3 / 1000;
// print the results to the textarea
print("\n z_1 = " + z_1);
print("\n z_2 = " + z_2);
print("\n z_3 = " + z_3);
print("\n z_4 = " + z_4);
print("\n z_5 = " + z_5);
]]></script>
<_-.XholonClass>
<!-- types of domain objects -->
<TheSystem/>
<SolarSystem/>
<Sun/>
<Planet>
<Earth/>
</Planet>
<Space/>
<TopOfAtmosphere/>
<Atmosphere/>
<Layer/> <!-- a vertical layer of the atmosphere, used for modeling -->
<Surface/>
<!-- quantities -->
<PhysicalProperty>
<Albedo/>
<!-- properties of the atmosphere -->
<Absorptivity/> <!-- a -->
<Emissivity/> <!-- e -->
<Temperature/> <!-- Kelvin -->
<Celsius/> <!-- Temperature in Celsius -->
<HeatCapacity/>
<SolarIn/> <!-- incoming solar radiation -->
<Ir/> <!-- a constant infrared parameter for the layer -->
<IrUp/> <!-- up-flowing infrared radiation -->
<IrDown/> <!-- down-flowing infrared radiation -->
<Z/> <!-- z "convective adjustment" -->
</PhysicalProperty>
<!-- constants -->
<Constant>
<StefanBoltzmannConstant/> <!-- Sigma -->
<SecondsPerTimeStep/>
<LapseMax/> <!-- "rate of temperature decrease with height" -->
<ZeroDegrees/> <!-- 0 degrees Celsius °C, expressed in degrees Kelvin K -->
</Constant>
<Constants/>
<Viewer/> <!-- this should be in the View part of the MVC structure -->
</_-.XholonClass>
<xholonClassDetails>
</xholonClassDetails>
<TheSystem>
<Constants>
<StefanBoltzmannConstant unit="W/m^2/K^4">5.6696e-8</StefanBoltzmannConstant>
<SecondsPerTimeStep unit="s">31536000.0</SecondsPerTimeStep> <!-- SecondsPerYear 365*24*60*60 -->
<LapseMax unit="K/km">6.5</LapseMax>
<ZeroDegrees unit="K">273.15</ZeroDegrees>
</Constants>
<SolarSystem>
<Sun/>
<Space/>
<Earth>
<TopOfAtmosphere roleName="tpf">
<Albedo>0.3</Albedo>
<SolarIn>0.0 W/m^2</SolarIn>
</TopOfAtmosphere>
<Atmosphere>
<!-- the first layer below is the bottommost layer of the Atmosphere -->
<Layer multiplicity="20">
<Temperature unit="K"> <!-- different initial temperature for each layer -->
<Celsius unit="°C"/>
</Temperature>
<Emissivity>0.34</Emissivity> <!-- "for infrared radiation" same value for each layer -->
<Absorptivity>0.0613</Absorptivity> <!-- "for solar radiation" same value for each layer -->
<HeatCapacity>1.0 J/m^3/K</HeatCapacity> <!-- not important in this model, so set it to 1 -->
<SolarIn>0.0 W/m^2</SolarIn>
<Ir>0.0 W/m^2</Ir>
<IrUp>0.0 W/m^2</IrUp>
<IrDown>0.0 W/m^2</IrDown>
<Z/>
</Layer>
</Atmosphere>
<Surface roleName="srf">
<Temperature unit="K">290.19
<Celsius unit="°C"/>
</Temperature>
</Surface>
</Earth>
</SolarSystem>
<Viewer/>
</TheSystem>
<Blockbehavior implName="lang:python:inline:"><![CDATA[
]]></Blockbehavior>
<Blockbehavior implName="lang:javascript:inline:"><![CDATA[
]]></Blockbehavior>
<TheSystembehavior implName="lang:webEditionjs:inline:"><![CDATA[
function postConfigure() {
$("textarea#btstrp_console").css("font-size", "12px");
$("textarea#btstrp_console").attr("rows", "3");
$("textarea#btstrp_console").attr("cols", "50");
}
]]></TheSystembehavior>
<Atmospherebehavior implName="lang:webEditionjs:inline:"><![CDATA[
function postConfigure() {
// initial temperatures from the base of the atmosphere up; T_1 to T_5; the temperature decreases as the altitude increases
//var initTempVal = [281.3, 267.0, 257.3, 245.6, 231.0]; // expected final values ?
var initTempVal = [273.15, 273.15, 273.15, 273.15, 273.15];
// z values; convective adjustments; z_1 to z_5
var zVal = [ 0.882, 2.9, 5.5, 9.2, 16.9];
// layers of the atmosphere from lowest to highest
var layers = this.parent().children(".Layer");
// get sigma which is a constant
var sigma = this.closest(".TheSystem").children(".Constants").children(".StefanBoltzmannConstant").attr("val");
var zeroDegrees = this.closest(".TheSystem").children(".Constants").children(".ZeroDegrees").attr("val");
// configure initial values of each layer
layers.each( function(index) {
$this = $(this);
var kelvin = $this.children(".Temperature");
var celsius = kelvin.children(".Celsius");
var ir = $this.children(".Ir");
var e = $this.children(".Emissivity").attr("val");
var z = $this.children(".Z");
kelvin.attr("val", initTempVal[index]);
celsius.attr("val", initTempVal[index] - zeroDegrees);
ir.attr("val", e * sigma * Math.pow(kelvin.attr("val"), 4));
z.attr("val", zVal[index]);
$this.attr("roleName", "lr" + (index+1));
});
}
]]></Atmospherebehavior>
<TopOfAtmospherebehavior implName="lang:webEditionjs:inline:"><![CDATA[
function postConfigure() {
// calculate incoming solar radiation
var thisTpf = this.parent();
var albedo = thisTpf.children(".Albedo").attr("val");
var solarIn = 1368 / 4 * (1 - albedo); // 239
thisTpf.children(".SolarIn").attr("val", solarIn);
}
]]></TopOfAtmospherebehavior>
<Layerbehavior implName="lang:webEditionjs:inline:"><![CDATA[
function act() {
var layerThis = this.parent();
var layerAbove = layerThis.next(".Layer");
if (!layerAbove.length) {
layerAbove = layerThis.closest(".Earth").children(".TopOfAtmosphere");
}
var layerBelow = layerThis.prev(".Layer");
if (!layerBelow.length) {
layerBelow = layerThis.closest(".Earth").children(".Surface");
}
//print("\n" + layerThis.attr("roleName") + " " + layerBelow.attr("roleName") + " " + layerAbove.attr("roleName") + " ");
var a = layerThis.children(".Absorptivity").attr("val");
var e = layerThis.children(".Emissivity").attr("val");
// inflows
var solarIn = layerAbove.children(".SolarIn").attr("val") * (1 - a);
var irFromAbove = layerAbove.children(".IrDown").attr("val") * (1 - e) + (layerThis.children(".Ir").attr("val")-0);
var irFromBelow = layerBelow.children(".IrUp").attr("val") * (1 - e) + (layerThis.children(".Ir").attr("val")-0);
layerThis.children(".SolarIn").attr("val", solarIn);
// temperature (a measure of the amount of energy in the layer)
// T_3(t) = T_3(t - dt) + (S_4_3 + IR_4_3 + IR_2_3 + C_2_3 - IR_3_4 - C_3_4 - S_3_2 - IR_3_2) * dt
var kelvinNode = layerThis.children(".Temperature");
var celsiusNode = kelvinNode.children(".Celsius");
var kelvin = kelvinNode.attr("val") + (solarIn + irFromAbove + irFromBelow) * dt;
// outflows
var irUp = irFromBelow * (1 - e) + (layerThis.children(".Ir").attr("val")-0);
var irDown = irFromAbove * (1 - e) + (layerThis.children(".Ir").attr("val")-0);
layerThis.children(".IrUp").attr("val", irUp);
layerThis.children(".IrDown").attr("val", irDown);
}
]]></Layerbehavior>
<Viewerbehavior implName="lang:webEditionjs:inline:"><![CDATA[
function postConfigure() {
this.xholoncreationservice('requireScript', 'raphael-min.js');
// create ports to the nodes that have the data for the viewer
var earth = this.closest(".TheSystem").children(".SolarSystem").children(".Earth");
this.data("tpf", earth.children(".TopOfAtmosphere"));
this.data("layers", earth.children(".Atmosphere").children(".Layer"));
this.data("srf", earth.children(".Surface"));
}
function act() {
var tpf = this.data("tpf");
var layers = this.data("layers");
var srf = this.data("srf");
if (this.application("getTimeStep") == 0) {
// use the existing SVG div, and remove any content it has
var mySVGDiv = $("div#mySVGDiv");
mySVGDiv.empty();
var x = 0, y = 0, width = 300, height = 350;
var shapeHeight = height / (layers.length + 2);
var linc = 0.7 / (layers.length + 2); // hsl color increment
var rgb = Raphael.getRGB("#283C5F"); // dark cornflowerblue
var hsl = Raphael.rgb2hsl(rgb.r, rgb.g, rgb.b);
var paper = Raphael(mySVGDiv[0], width, height);
var hex = Raphael.hsl(hsl.h, hsl.s, hsl.l);
// top of atmosphere (tpf)
var vtpf = paper.rect(x, y, width, shapeHeight).attr("fill", hex).attr("stroke", hex);
vtpf.mouseover( function() {
print("\n" + tpf.attr("roleName"));
});
this.data("vtpf", vtpf);
// atmosphere layers
layers.each( function(index) {
y += shapeHeight;
hsl.l = hsl.l + linc;
hex = Raphael.hsl(hsl.h, hsl.s, hsl.l);
var shape = paper.rect(x, y, width, shapeHeight).attr("fill", hex).attr("stroke", hex);
var thisLayer = $(this);
shape.mouseover( function() {
print("\n" + thisLayer.attr("roleName") + " " + thisLayer.children(".Temperature").attr("val"));
});
// inter-layer flux shapes
var yTop = y - 5;
var yBot = y + 5;
paper.path("M 25 " + yTop + " L 35 " + yTop + " L 30 " + yBot + " Z").attr("fill", "yellow").attr("stroke", "yellow"); // solarIn
if (index > 0) {
paper.path("M 50 " + yTop + " L 60 " + yTop + " L 55 " + yBot + " Z").attr("fill", "red").attr("stroke", "red"); // irDown
}
paper.path("M 75 " + yBot + " L 85 " + yBot + " L 80 " + yTop + " Z").attr("fill", "red").attr("stroke", "red"); // irUp
paper.path("M 100 " + yBot + " L 110 " + yBot + " L 105 " + yTop + " Z").attr("fill", "green").attr("stroke", "green"); // convection
});
y += shapeHeight;
hsl.l = hsl.l + linc;
hex = Raphael.hsl(hsl.h, hsl.s, hsl.l);
// surface (srf)
var vsrf = paper.rect(x, y, width, shapeHeight).attr("fill", hex).attr("stroke", hex);
vsrf.mouseover( function() {
print("\n" + srf.attr("roleName") + " " + srf.children(".Temperature").attr("val"));
});
this.data("vsrf", vsrf);
// inter-layer flux shapes
var yTop = y - 5;
var yBot = y + 5;
paper.path("M 25 " + yTop + " L 35 " + yTop + " L 30 " + yBot + " Z").attr("fill", "yellow").attr("stroke", "yellow"); // solarIn
paper.path("M 50 " + yTop + " L 60 " + yTop + " L 55 " + yBot + " Z").attr("fill", "red").attr("stroke", "red"); // irDown
paper.path("M 75 " + yBot + " L 85 " + yBot + " L 80 " + yTop + " Z").attr("fill", "red").attr("stroke", "red"); // irUp
paper.path("M 100 " + yBot + " L 110 " + yBot + " L 105 " + yTop + " Z").attr("fill", "green").attr("stroke", "green"); // convection
}
}
]]></Viewerbehavior>
<Blockbehavior implName="lang:bsh:inline:"><![CDATA[
]]></Blockbehavior>
<Blockbehavior implName="lang:jruby:inline:"><![CDATA[
]]></Blockbehavior>
<Blockbehavior implName="lang:groovy:inline:"><![CDATA[
]]></Blockbehavior>
<SvgClient><Attribute_String roleName="svgUri"><![CDATA[data:image/svg+xml,
<svg width="100" height="50" xmlns="http://www.w3.org/2000/svg">
<g>
<title>Radiative-Convective Model</title>
<rect id="TheSystem/Atmosphere" fill="skyblue" height="50" width="50" x="25" y="0"/>
</g>
</svg>
]]></Attribute_String><Attribute_String roleName="setup">${MODELNAME_DEFAULT},${SVGURI_DEFAULT}</Attribute_String></SvgClient>
</XholonWorkbook>