# Circulating Cells, an atmospheric convection simulation.
#
# Copyright (C) 2011-2012 Kevan Hashemi, kevan@hashemifamily.com
#
set p(name) "Circulating Cells, Version 11.2"
#
# This program simulates the transport of heat in an array of atmospheric gas
# cells, with surface blocks of either water or sand. Water evaporates from
# the surface blocks, condenses in into clouds in cool gas cells, turns into
# snow in freezing gas cells, descends as snowflakes, melts into drops
# in warm cells, and falls to the surface as rain. The surface blocks are
# warmed by the light of the sun. The surface gas warms in contact with the
# surface. The warm cells rise. As they rise, their pressure drops. They
# expand. As they expand they cool. Nevertheless, they continue to rise because
# they are warmer than their neighbors. If water vapor condenses within the
# cell this further warms them so that they rise faster. When they reach the
# top of the array, they radiate heat into space. We call this top row of the
# array the "tropopause" because it is the point in our simulated atmosphere
# where convection stops. The gas of the asmosphere is partly transparent to
# the long-wave radiation emitted by the surface blocks. Clouds, on the other
# hand, are opaque to long-wave radiation. A cloud radiates like a black
# body. A gas cell radiates like a gray body.
#
# To run the program, start by downloading this file and saving to disk. On Linux
# or MacOS open a terminal shell, navigate to the same directory as the program,
# and type "wish CC_n.txt" where n is the version of this program you want to run.
# Version 1 is called CC_1.tcl. On Windows, you have to download the program that
# will take this file and turn it into the simulation process. Go to the Active
# State website and download ActiveTcl for Windows. It's free, and you don't have
# to sign anything or fill any forms if you don't want to when you download.
# Install the program and run it. In the File menu, select Source. A file browser
# will open. Select the CC_n.tcl file. The simulation will start.
#
# The program divides the column into cells each containing the same mass of gas.
# The column of air is so high that the pressure at the bottom significantly
# greater than at the top. Air rising expands and cools. Air that descends to make
# way for rising air will compress and heat.
#
# The comments in the code below explain how we calculate the cell heights for the
# display, the center of mass of cells for the circulation, and the temperature
# and pressure of the cells as we swap them around. The basis of all these
# calculations is the equal mass of the cells, which means that the drop in
# pressure going up one cell height is constant in absolute terms, not in relative
# terms. If the cell mass is 100 kg/m2, this is a weight of 1000 N/m2 so the pressure
# just below a cell will be 1 kPa greater than just above. We assume that the
# cells are small enough that we can use the pressure at the bottom as the average
# pressure of the entire cell.
#
# The simulation of convection is random in that it tests one randomly-chosen block
# of cells after another to see if circulation will generate enough work to cause
# the cells to rotate. The cells are color-coded according to temperature. Blue is
# cool and red is hot, as shown in the legend at the top. You can reset the cells
# with the Reset button. As the program runs, it displays the number of random
# blocks it has worked on since the run began. We refer to each choice of block as
# an iteration of the program. In an array of 10x10 cells, each cell is part of
# a block every 25 iterations on average.
#
# Convection moves cells around. You can mark any cell you like by clicking on it.
# The border of the cell will turn from white to black and you can watch it move
# around. You can mark as many cells as you like. Click a marked cell a second
# time and it will return to normal. You can obtain cell data in text form with
# the report function and the Data Window. The default report function prints
# out the time in hours and the average temperature of each row.
#
# Cells can contain water vapor, water droplets, and snow flakes. Water vapor
# enters a cell in one of two ways. It can evaporate from a water surface block or
# it can evaporate from a droplet when the temperature is above freezing. Droplets
# form in only one way: by condensation of water vapor when the temperature of the
# cell is above the freezing point of water. Droplets vanish in one of three ways.
# They can evaporate into water vapor when the temperature is above freezing, or
# they can be transformed into snow flakes by the Bergeron process when the
# temperature is below freezing. In the Bergeron process, snow flakes start to
# form in the gas, which is saturated with water vapor on account of its
# containing water droplets. The snow flakes remove the water vapor and the
# droplets evaporate. Our simulation skips the intermediate water vapor stage, and
# simply has the water droplets turning into snow flakes. The snow flakes disappear
# from the cell by falling out of it. If they enter a cell where the temperature is
# above freezing, they will melt. Once melted, they are rain drops, and these fall
# into the cell below after cooling down the surrounding gas cell by taking from it
# their latent heat of fusion. The speed of snow and rain fall are set in meters
# per second by snow_speed_mps and rain_speed_mps. The speed with which droplets
# turn into snow is freeze_rate_gps in grams of water per kilogram of dry air
# per second. The speed with which snow flakes melt is melt_rate_gps. The third
# way water droplets can vanish is by sinking, which they do at cloud_speed_mps.
#
# Any cell with more than wc_opaque in g/kg of liquid or solid water is opaque to
# to long-wave radiation. Such a cell radiates like a black body. We ignore the
# absorption of long-wave radiation by water vapor. Cells with less liquid and
# solid water we treat as if they were dry. Our gas is transparent to a fraction
# of the long-wave spectrum, as specified by the transparency_fraction. At other
# wavelengths it is opaque. There is no wavelength at which it is somewhere
# in-between. Thus the radiatioin emitted by a gas cell is absorped entirely by
# the gas cells above and below. The surface blocks we assume are opaque to both
# short-wave radiation from the sun and all long-wave radiation. A dry gas cell can
# radiate heat only when it arrives at the tropopause (the top row). But a wet gas cell
# can radiate heat at any altitude, because it radiates at wavelengths the gas above
# and below do not absorb.
#
# This program is free software; you can redistribute it and/or modify it under
# the terms of the GNU General Public License as published by the Free Software
# Foundation; either version 2 of the License, or (at your option) any later
# version.
#
# Define the dimensions of the cell array.
set p(array_width) 30
set p(array_height) 15
# The initial temperature of all cells in Kelvin.
set p(T_initial) 280
# The pressure of the top and bottom cells in the column in Pascal.
set p(p_bottom_kPa) 100.0
set p(p_top_kPa) 50.0
# The pressure step from one cell row to the next is constant, which means that
# the mass of the cells is equal, even as we rise up the column.
set p(p_step_kPa) [expr 1.0 * ($p(p_bottom_kPa) - $p(p_top_kPa)) / $p(array_height)]
set p(p_step_fraction) [expr 1.0 * $p(p_step_kPa) / $p(p_bottom_kPa)]
# When air expands adiabatically, its pressure and temperature are related by
# p(1-Cp/Cv)*T^(Cp/Cv)=constant, where Cp and Cv are the heat capacities of air
# at constant pressure and constant volume respectively. We also have Cp=Cv+R
# where R is the constant pV/T for 1 kg of air. Thus we have the relation
# T2/T1 = (p2/p1)^(R/Cp). Units for R are J/K and for Cp are J/kgK.
set p(R_gas) 287
set p(Cp_gas) 1003
# And for water.
set p(R_water) 462
set p(Cp_water) 1840
# We choose a value for gravitational acceleration throughout the array, which
# we express in units of N/kg, or m/s/s.
set p(g) 10.0
# The coefficient of convective heat transfer from a surface to the gas above,
# in Watt per square meter per degree centigrade. The coefficient applies only
# when the surface is warmer than the gas.
set p(convection_coeff) 20
# Here is a routine to calculate the saturation concentration of water vapor
# in terms of the absolute pressure of the air in which it is dissolved. The
# temperature is in Kelvin. We specify the pressure with the row number of the
# gas cell, which allows us to calculate the pressure in the routine, and thus
# simplify the calling routine. The result is saturation concentration in grams
# of water vapor dissolved in each kilogram of dry air. We use an approximation
# we derived in our Condensation Point post.
proc saturation_concentration {Tc j} {
global p
if {$Tc>250.0} {
set pc [expr $p(p_bottom_kPa) - $j * $p(p_step_kPa)]
return [expr ($Tc-250.0)*($Tc-250.0)*$pc/8000.0]
} {
return 0.0
}
}
# Suppose we have gas with water vapor concentration X. If X is greater than
# the saturation concentration, vapor must condense. But as it condenses, the
# vapor releases its latent heat of evaporation, and so warms uup the gas, which
# increases its saturation concentration. Conversely, suppose X is less than the
# saturation concentration, and we have water droplets in the air. The droplets
# must evaporate, but as they do, they remove their latent heat of eveporation
# from the gas, and so cool the gas down and decrease the saturation concentration.
# The following routine determines the temperature at which the gas will, by either
# process, be saturated with water vapor. This temperature is the dew point, which
# we refer to at Ts, the saturation temperature, in the routine and other places.
# The combination of our saturation concentration equation and the linear heating
# relation gives us a quadratic equation that we solve with the quadratic formula,
# choosing the positive root if it exists. For temperatures below 250K, where
# our approximation switches to zero saturation concentration, we solve
# for the balance point more simply.
proc dew_point {Tc Xc j} {
global p
set pc [expr $p(p_bottom_kPa) - $j * $p(p_step_kPa)]
set b [expr 8000.0*$p(Cp_gas)/$pc/$p(Qe_water)-2*250.0]
set c [expr 250.0*250.0-8000.0*$p(Cp_gas)/$pc/$p(Qe_water)*$Tc-8000.0*$Xc/$pc]
if {$b*$b-4*$c>0} {
set Ts [expr (-$b+sqrt($b*$b-4*$c))/2.0]
if {$Ts<250.0} {
set Ts [expr $Tc + $Xc*$p(Qe_water)/$p(Cp_gas)]
}
} {
set Ts [expr $Tc + $Xc*$p(Qe_water)/$p(Cp_gas)]
}
return $Ts
}
# Here is a routine to calculate the evaporation rate of water from the surface
# of a lake or ocean. The temperature is that of the water in Kelvin. The
# pre-existing concentration of water vapor in the air is in g/kg. The evaporation
# rate is in gram per square meter per second. We use an approximation we derived
# in our Evaporation Rate post.
proc evaporation_rate {T x} {
global p
return [expr (($T-250.0)*($T-250.0)/80.0-$x)/40.0]
}
# Latent heat of evaporation and fusion for water, in Joules per gram. Freezing
# point of water droplets in Kelvin.
set p(Qe_water) 2200
set p(Qf_water) 330
# The freezing point of cloud droplets is the temperature at which snow flakes
# start to form, absorbing existing water vapor and leading to the evaporation
# of cloud droplets to maintain water vapor at its saturation level. We set this
# temperature a few degrees below the freezing point of water, because spontaneous
# crystal formation does not occur at exactly 0C.
set p(Tf_droplets) 268
# The melting point of snow flakes is sligthly higher than the melting point of ice
# because the melting takes time and also cools the surrounding air by absorption
# of latent heat of fusion.
set p(Tm_ice) 278
# Here we specify the reflectivity of clouds, in terms of cloud thickness. We
# measure cloud thickness as the depth of the water layer that would be made
# by its combined droplets. We specify the penetration depth of sunlight in
# a cloud, which is the cloud depth in millimeters that will reflect 63% of
# sunlight.
set p(Lc_water) 3.0
# The minimum water droplet concentration in g/kg that will be be marked as a
# cloud. The mass of cloud water per square meter is this concentration multipled
# by the cell mass. For 300-kg cells, a concentration of 1.0 g/kg gives 300 g/m2,
# which is 0.3 mm depth, which reflects roughly 10% of incoming sunlight when
# Lc_water is 3.0 mm.
set p(wd_cloud) 1.0
# The minimum rain drop concentration in g/kg that will be marked as a rain shower,
# provided the cell is not a cloud. The depth of the rain we calculate in the same
# way as for clouds. It the cells are 300 kg/m2 and the rain drop concentraion is
# 0.3 g/kg, the cell contains 30 g/m2 or 0.09 mm depth of water. If this rain
# falls at 5 m/s, it will all drop out of the cell in roughly 100 s, so we have
# 3 mm/hr of rainwater, which is light rain.
set p(rd_shower) 0.3
# The minimum snow flake concentration in g/kg that will be marked as a snow
# fall, provided it is not a cloud nor a rain shower.
set p(sf_snowfall) 0.1
# The threshold concentration in g/kg of liquid and solid water at which we
# consider a cell to be opaque to long-wave radiation. A concentration of
# 0.3 g/kg in a 300-kg/m2 cell represents a depth of 90 um, which is sufficient
# to absorb 98% of all long-wave radiation.
set p(wc_opaque) 0.3
# Ice crystals float down, rain falls down, and water droplets sink very slowly.
# We specify each speed here.
set p(cloud_speed_mps) 0.003
set p(snow_speed_mps) 1.0
set p(rain_speed_mps) 5.0
# If the temperature of a cloud is below Tf, snow starts to form. We assume a
# rate in grams per kilogram per second. For snow melting above Tm, we have another
# rate in grams per kilogram per second.
set p(freeze_rate_gps) 0.001
set p(melt_rate_gps) 0.01
# Heat capacity of various surface blocks, in degree centigrade per Joule per
# square meter.
set p(heat_capacity_sand) 140000.0 ;# Ten-centimeter layer of sand
set p(heat_capacity_water) 4200000.0 ;# One meter layer of water
# Simulation timing.
set p(iterations_per_s) 1.0
set p(s_per_hr) 3600
set p(day_length_hr) 24
set p(mid_day_hr) 12.0
set p(mid_night_hr) 0.0
set p(daylight_fraction) 0.5
set p(solar_time) 0.0
# Q_sun is in W/m2 for the surface material cells during the daylight hours.
# The Q_sun parameter contains the current heating rate of the sun, for the purpose
# of display or reporating.
set p(Q_sun) 350.0
set p(Q_current) 0.0
# These two variables will hole the upwelling and downwelling radiation average
# values in W/m2. The upwelling radiation is the total long-wave radiation power
# leaving the atmosphere in the upward direction, and passing into space. It
# includes the radiation from the planet surface. The downwelling radiation is
# the total long-wave radiation power received by the planet surface from the
# atmosphere.
set p(Q_upwelling) 0.0
set p(Q_downwelling) 0.0
# We use stefan's constant to calculate radiation from various surfaces. It's
# units are W/m2/K4. We use the transparency fraction to determine the fraction
# of long-wave radiation radiated by a black body that will pass through the
# atmospheric gas without clouds. We use one minus this fraction to determine the
# heat radiated by the tropopause, which is the top layer of cells, and is itself
# transparent to the same wavelenthgs, but unlike the lower layers of the atmosphere
# the tropopause has no further opaque atmosphere above to stop its radiation propagating
# out into space.
set p(stefan_const) 5.7e-8
set p(transparency_fraction) 0.5
# Here we set the depth of liquid water and ice in a cell that causes it to be
# opaque to all long-wave radiation. Cells with less water we allow to radiate
# like dry cells.
set p(opaque_depth_mm) 0.1
# Control variables.
set p(control) "Stop"
set p(time_control) "Day"
set p(accelerate) 1
set p(tracking) 0
set p(reporting) 1
set p(data_description) ""
set p(counter) 0
set p(report_interval) 3600
set p(plot_interval) 1000
set p(update_interval) 100
set p(check_interval) 10000
set p(heating_interval) 100
# Constants for display of temperatures.
set p(max_T) 320
set p(min_T) 240
# Constants that control behavior of marked cells.
set p(mark_num) 0
# Math constants.
set p(pi) 3.141592654
# We calculate the cell mass per square meter of its bottom surface area. The
# approximate cell height is the mass divided by the density at the initial
# temperature at the bottom of the array.
set p(cell_mass_kg) [format %.1f [expr 1000.0 * $p(p_step_kPa) / $p(g)]]
set p(cell_height_m) [format %.1f \
[expr $p(cell_mass_kg) * $p(T_initial) * $p(R_gas) / $p(p_bottom_kPa) / 1000.0]]
# The program's selection interval is the average time between selections of
# a cell for circulation and water balancing. Four cells are selected on each
# iteration of the program.
set p(selection_interval_s) [expr $p(array_width) * $p(array_height) \
/ $p(iterations_per_s) / 4.0]
# The impetus threshold is the minimum impetus for circulation at which we will
# allow a rotation of four cells to take place. The threshold is a function of
# the cell height and the array area. The time between oppertunities for a hot
# cell to rise is roughly equal to twice the selection interval. If we assume
# the cells accelerate at a constant rate to their maximum speed as a result
# of the circulation impetus, then decelerate at a constant rate by the action
# of viscocity, their average speed will be half the maximum speed. When they
# are moving at their maximum speed, their kinetic energy can be no greater
# than the impetus for circulation, and half the maximum speed must be adequate
# to travel one cell height in the available time. We calculate a likely value
# for the impetus threshold at start-up, but the user can change this value
# later in the configuration array. We set the threshold at twice the minimum
# energy for circulation.
set p(impetus_threshold) [format %.1f [expr 2.0 * 0.5 * \
pow(0.5 * $p(cell_height_m) / $p(selection_interval_s), 2.0)]]
# Calculate display dimensions in units of pixels.
set p(rectangle_width) 25
set p(rectangle_height) [expr round(1.2*$p(rectangle_width))]
set p(canvas_width) [expr $p(rectangle_width) * $p(array_width)]
set p(canvas_height) [expr $p(rectangle_height) * ($p(array_height)+1)]
# A list of parameters we want displayed in the configuration panel and saved to
# our array test files.
set p(config_names) "Q_sun transparency_fraction day_length_hr daylight_fraction \
iterations_per_s T_initial impetus_threshold \
report_interval plot_interval update_interval heating_interval \
counter Cp_gas R_gas heat_capacity_sand heat_capacity_water \
Lc_water Tf_droplets Tm_ice convection_coeff wd_cloud rd_shower \
sf_snowfall snow_speed_mps rain_speed_mps cloud_speed_mps \
freeze_rate_gps melt_rate_gps"
# Configure the windows.
set p(main_window) .
set p(data_window) .data
wm title . $p(name)
# The following three routines turn the standard input into a TclTk console. At
# the console, we can enter Tcl commands to manipulate and query the simulation.
# We use these routines only if the TclTk we're using does not provide its
# own console with the "console" command.
proc console_start {} {
fileevent stdin readable console_execute
fconfigure stdin -buffering line
console_prompt
}
proc console_prompt {} {
puts -nonewline stdout "CC% "
flush stdout
}
proc console_execute {} {
gets stdin line
catch {uplevel $line} result
if {$result != ""} {puts stdout $result}
console_prompt
}
# Create a frame for buttons at the top.
frame .f
pack .f -side top -fill x
foreach a {Control Counter} {
set b [string tolower $a]
label .f.$b -textvariable p($b) -width 8
pack .f.$b -side left
}
foreach a {Run Stop Reset} {
set b [string tolower $a]
button .f.$b -text $a -command [list set p(control) $a]
pack .f.$b -side left -expand 1
}
foreach a {Configure Data Save Load} {
set b [string tolower $a]
button .f.$b -text $a -command $b
pack .f.$b -side left -expand 1
}
# We create a console button if the TclTk we're using provides its
# own console. Otherwise we turn the standard input into a terminal.
if {[info command console] != ""} {
button .f.console -text Console -command "console show"
pack .f.console -side left -expand 1
} {
console_start
}
frame .c
pack .c -side top -fill x
label .c.tl -text "Phase (hr)" -fg orange
label .c.tv -textvariable p(solar_time) -width 4 -fg orange
pack .c.tl .c.tv -side left
label .c.hl -text "Sun (W/m2)" -fg orange
label .c.hv -textvariable p(Q_current) -width 5 -fg orange
pack .c.hl .c.hv -side left
tk_optionMenu .c.config p(time_control) "Day" "Night" "Cycle"
pack .c.config -side left -expand 1
foreach a {Accelerate Tracking Reporting} {
set b [string tolower $a]
checkbutton .c.$b -variable p($b) -text $a
pack .c.$b -side left -expand 1
}
# Create the legend.
canvas .l -height 20 -width $p(canvas_width)
pack .l -side top -fill x
# This procedure draws the temperature legend on the main window,
# defines the relationship between temperature and color in the cell display.
proc temperature_scale {} {
global p
# Define the blue, green, and red temperatures in terms of the
# maximum and minimum temperatures.
set p(blue_T) $p(min_T)
set p(green_T) [expr round(0.5*($p(max_T)+$p(min_T)))]
set p(red_T) $p(max_T)
# Define or re-define a routine to convert temperature into color.
set p(colorreach) [expr ($p(red_T) - $p(blue_T)) / 2.0]
proc get_temp_color {T} {
global p
if {$T < $p(blue_T)} {set T $p(blue_T)}
if {$T > $p(red_T)} {set T $p(red_T)}
foreach rgb "red green blue" {
set $rgb [expr round(255.0*exp(\
-pow(1.0*([set p($rgb\_T)]-$T)/$p(colorreach),2)))]
}
return [format "#%02x%02x%02x" $red $green $blue]
}
# Draw the legend with the help of our new conversion routine.
for {set T $p(blue_T)} {$T <= $p(red_T)} {set T [expr $T+1]} {
set x [expr $p(canvas_width) * (0.1 + 0.8 * ($T-$p(blue_T)) / ($p(red_T)-$p(blue_T)))]
set y [expr 5]
.l create rectangle [expr $x-4] [expr $y-2] [expr $x+4] [expr $y+2] \
-outline [get_temp_color $T] -fill [get_temp_color $T]
if {[expr fmod($T,10)] == 0} {
.l create text $x [expr $y + 10] -text [expr round($T)]
}
}
}
# Draw the temperature scale.
temperature_scale
# Create a place for the simulation display.
canvas .d -height [expr $p(canvas_height)] -width [expr $p(canvas_width)]
pack .d -side top
# Create the empty cell array that we will shortly fill up.
# Create an array of cells and an accompanying array of rectangles in the
# display. Each cell has a rectangle and a unique array index (j,i) where
# j is the row and i is the column. Row 0 is the bottom row. Column 0 is
# the leftmost column. The array index, the temperature of the cell, the
# name of the display rectangle, and the type of display for the cell are
# all given by the cell element in the cell list. In the display, we
# give the cells in each row a height that is in proportion to the volume
# occupied by the cell masses at that altitude. If the average cell height
# is H, the bottom cell height is H-b and the top is H+b, where the ratio
# (H+b)/(H-b) = p_bottom/p_top.
set cells ""
set y_bottom $p(rectangle_height)
set a [expr 1.0 * $p(p_bottom_kPa) / $p(p_top_kPa)]
set b [expr $p(rectangle_height) * ($a - 1.0) / ($a + 1.0)]
for {set j 0} {$j < $p(array_height)} {incr j} {
set y_top [expr $y_bottom + ($p(rectangle_height)-$b)+2.0*$b*$j/$p(array_height)]
set row ""
for {set i 0} {$i < $p(array_width)} {incr i} {
set left [expr $i*$p(rectangle_width)]
set top [expr $p(canvas_height) - $y_top + 1]
set right [expr ($i+1)*$p(rectangle_width) - 1]
set bottom [expr $p(canvas_height) - $y_bottom]
set name [.d create rectangle \
$left $top $right $bottom \
-outline white -fill [get_temp_color $p(T_initial)]]
.d bind $name {mark_cell %x %y %W}
lappend row "$j $i $name unmarked $p(T_initial) 0.0 0.0 0.0 0.0"
}
lappend cells $row
set y_bottom $y_top
}
# We create a row of surface blocks to represent the surface of the planet.
# The i'th block is beneath gas cell i in the bottom row of gas cells.
set blocks ""
set y_bottom 0
set y_top [expr $p(rectangle_height)]
for {set i 0} {$i < $p(array_width)} {incr i} {
set left [expr $i*$p(rectangle_width)]
set top [expr $p(canvas_height) - $y_top + 1]
set right [expr ($i+1)*$p(rectangle_width) - 1]
set bottom [expr $p(canvas_height) - $y_bottom]
set name [.d create rectangle \
$left $top $right $bottom \
-outline white -fill [get_temp_color $p(T_initial)]]
.d bind $name {configure_block %x %y %W}
lappend blocks "$i $name sand $p(T_initial)"
}
# Routines to get and set the temperature of a cell or block.
proc get_temp {j i} {
global cells
return [format %.3f [lindex $cells $j $i 4]]
}
proc set_temp {j i T} {
global cells
lset cells $j $i 4 [format %.3f $T]
return $T
}
proc get_block_temp {i} {
global blocks
return [format %.3f [lindex $blocks $i 3]]
}
proc set_block_temp {i T} {
global blocks
lset blocks $i 3 [format %.3f $T]
return $T
}
# Routine to calculate average temperature of the cells in a row.
proc get_row_temp {j} {
global p
set ave 0.0
for {set i 0} {$i < $p(array_width)} {incr i} {
if {$j < 0} {
set T [get_block_temp $i]
} {
set T [get_temp $j $i]
}
set ave [expr $ave + $T]
}
return [format %.1f [expr 1.0 * $ave / $p(array_width)]]
}
# Routines to get and set the type of a cell or block.
proc get_type {j i} {
global cells
return [lindex $cells $j $i 3]
}
proc set_type {j i type} {
global cells
lset cells $j $i 3 $type
return $type
}
proc get_block_type {i} {
global blocks
return [lindex $blocks $i 2]
}
proc set_block_type {i type} {
global blocks
lset blocks $i 2 $type
return $type
}
# Routines to get and set the water vapor concentration of a cell, in units
# of grams of water per kilogram of dry air.
proc get_wvc {j i} {
global cells
return [format %.3f [lindex $cells $j $i 5]]
}
proc set_wvc {j i c} {
global cells
lset cells $j $i 5 $c
return $c
}
# Routines to get and set the cloud droplet concentration of a cell.
proc get_cdc {j i} {
global cells
return [format %.3f [lindex $cells $j $i 6]]
}
proc set_cdc {j i c} {
global cells
lset cells $j $i 6 $c
return $c
}
# Routines to get and set the snow flake concentration of a cell.
proc get_sfc {j i} {
global cells
return [format %.3f [lindex $cells $j $i 7]]
}
proc set_sfc {j i c} {
global cells
lset cells $j $i 7 $c
return $c
}
# Routines to get and set the rain drop concentration of a cell.
proc get_rdc {j i} {
global cells
return [format %.3f [lindex $cells $j $i 8]]
}
proc set_rdc {j i c} {
global cells
lset cells $j $i 8 $c
return $c
}
# Routine to get the name of a cell or block.
proc get_name {j i} {
global cells
return [lindex $cells $j $i 2]
}
proc get_block_name {i} {
global blocks
return [lindex $blocks $i 1]
}
# Routine to calculate the altitude of a row, which we take
# to be half-way up the height of that row, plus the full height
# of all the rows underneath. We use the row temperatures from
# the current cell array, and the fixed row pressures to determine
# the average height of each row. The altitude is in meters.
proc get_row_altitude {j} {
global p
set altitude 0
set pressure $p(p_bottom_kPa)
for {set jj 0} {$jj <= $j} {incr jj} {
set h [expr 1.0 * $p(cell_mass_kg) * $p(R_gas) * [get_row_temp $j] \
/ $pressure / 1000]
set pressure [expr $pressure - $p(p_step_kPa)]
set altitude [expr $altitude + $h]
}
return [format %.1f [expr $altitude - $h/2.0]]
}
# water_balance looks at a cell and allows water droplets to evaporate,
# water vapor to condense, snow flakes for form, snow flakes to
# melt, and rain drops to fall. In our variables, the suffix "c"
# indicates the gas cell, and "s" the saturation state. We note that
# the units for Qe_water, the latent heat of evaporation for water,
# are in Joules per gram of water, while the water droplet and water
# vapor concentrations are in grams of water per kilogram of dry gas.
# The specific heat capacity of the gas at constant pressure is Cp_gas,
# which is in units of Joules per kilogram of dry gas. We ignore the
# contribution to heat capacity made by the water vapor and water
# droplets, their concentration being less than 1% most of the time.
proc water_balance {j i} {
global p
# Extract the characteristics of the cell into intermediate variables
# for our convenience.
set Tc [get_temp $j $i]
set Xc [get_wvc $j $i]
set Dc [get_cdc $j $i]
set Fc [get_sfc $j $i]
set Rc [get_rdc $j $i]
# If cloud droplets exist, allow some of them to sink into the cell below
# or land upon the surface.
if {$Dc > 0} {
set sinking_fraction [expr \
$p(cloud_speed_mps) * $p(selection_interval_s) / $p(cell_height_m)]
if {$sinking_fraction > 1.0} {set sinking_fraction 1.0}
set loss [expr $Dc*$sinking_fraction]
set Dc [expr $Dc-$loss]
if {$j>0} {
set_cdc [expr $j-1] $i [expr [get_cdc [expr $j-1] $i] + $loss]
}
}
# If rain drops exist, let some of it fall to the cell below or land
# upon the surface.
if {$Rc > 0} {
set sinking_fraction [expr \
$p(rain_speed_mps) * $p(selection_interval_s) / $p(cell_height_m)]
if {$sinking_fraction > 1.0} {set sinking_fraction 1.0}
set loss [expr $Rc*$sinking_fraction]
set Rc [expr $Rc-$loss]
if {$j>0} {
set_rdc [expr $j-1] $i [expr [get_rdc [expr $j-1] $i]+$loss]
}
}
# If snow flakes exist, and the cell is warm enough, we allow them
# to melt and turn into rain drops. They don't fall to the surface
# until the next time this cell is balanced.
if {($Tc >= $p(Tm_ice)) && ($Fc > 0)} {
set melt [expr $p(melt_rate_gps) * $p(selection_interval_s)]
if {$melt > $Fc} {set melt $Fc}
set Tc [expr $Tc - $melt*$p(Qf_water)/$p(Cp_gas)]
set Rc [expr $Rc + $melt]
set Fc [expr $Fc - $melt]
}
# If snow flakes still exist, allow them to fall into the cell below.
# This is the cloud snowing upon the cell below. If the cell below is a
# surface block, take the latent heat of fusion of the ice from the surface
# block.
if {$Fc>0.0} {
set sinking_fraction [expr \
$p(snow_speed_mps) * $p(selection_interval_s) / $p(cell_height_m)]
if {$sinking_fraction > 1.0} {set sinking_fraction 1.0}
set loss [expr $Fc*$sinking_fraction]
set Fc [expr $Fc-$loss]
if {$j>0} {
set_sfc [expr $j-1] $i [expr [get_sfc [expr $j-1] $i] + $loss]
} {
set Tb [get_block_temp $i]
set Cb [set p(heat_capacity_[get_block_type $i])]
set_block_temp $i [expr $Tb - $p(cell_mass_kg) * $loss * $p(Qf_water) / $Cb]
}
}
# We establish balance between water vapor, water droplets
# and gas temperature, ignoring any tendancy of the water vapor
# to deposit into snow flakes.
set Ts [dew_point $Tc $Xc $j]
set Xs [saturation_concentration $Ts $j]
if {$Xc > $Xs} {
set Tc $Ts
set Dc [expr $Dc + ($Xc-$Xs)]
set Xc $Xs
} elseif {$Dc > 0.0} {
if {$Dc <= ($Xs-$Xc)} {
set Tc [expr $Tc - $Dc*$p(Qe_water)/$p(Cp_gas)]
set Xc [expr $Xc + $Dc]
set Dc 0.0
} else {
set Tc [expr $Tc - ($Xs-$Xc)*$p(Qe_water)/$p(Cp_gas)]
set Dc [expr $Dc - ($Xs-$Xc)]
set Xc $Xs
}
}
# If the cell is cold enough, we assume snow flakes form rapidly
# and absorb water vapor, which in turn leads to the evaporation
# of all cloud droplets. The net result is that the a fraction
# of the cloud droplets turn into ice and give up their latent heat
# of fusion. This is the formation of snow. The snow flakes don't
# fall into the cell below until the next time this cell is balanced.
if {$Tc < $p(Tf_droplets)} {
set freeze [expr $p(freeze_rate_gps) * $p(selection_interval_s)]
if {$freeze > $Dc} {set freeze $Dc}
set Tc [expr $Tc + $freeze*$p(Qf_water)/$p(Cp_gas)]
set Fc [expr $Fc + $freeze]
set Dc [expr $Dc - $freeze]
}
# Assert the new cell charactersistics.
set_cdc $j $i $Dc
set_wvc $j $i $Xc
set_temp $j $i $Tc
set_sfc $j $i $Fc
set_rdc $j $i $Rc
}
# swap exchanges the marking and physical properties of cell
# j i with cell jj ii, leaving the display name and coordinates
# of the cell intact in the array list element that represents
# the cell.
proc swap {j i jj ii} {
global cells
set saved [lrange [lindex $cells $j $i] 3 end]
lset cells $j $i "$j $i [get_name $j $i] [lrange [lindex $cells $jj $ii] 3 end]"
lset cells $jj $ii "$jj $ii [get_name $jj $ii] $saved"
}
# To plot a cell, we fill its rectangle with a solid color that indicates
# its temperature, or in the case of a surface block we use the outline to
# indicate the material.
proc plot {j i} {
global p
set Tc [get_temp $j $i]
set name [get_name $j $i]
set type [get_type $j $i]
set Fc [get_sfc $j $i]
set Dc [get_cdc $j $i]
set Rc [get_rdc $j $i]
switch -glob $type {
"unmarked" {set outline white}
"mark*" {set outline black}
}
if {$Dc >= $p(wd_cloud)} {
.d itemconfigure $name -fill gray50 -outline $outline
} elseif {$Rc > $p(rd_shower)} {
.d itemconfigure $name -fill gray25 -outline $outline
} elseif {$Fc > $p(sf_snowfall)} {
.d itemconfigure $name -fill gray100 -outline $outline
} else {
.d itemconfigure $name -fill [get_temp_color $Tc] -outline $outline
}
}
proc plot_block {i} {
set Tb [get_block_temp $i]
set name [get_block_name $i]
set type [get_block_type $i]
switch -glob $type {
"water" {.d itemconfigure $name -fill [get_temp_color $Tb] -outline blue}
"sand" {.d itemconfigure $name -fill [get_temp_color $Tb] -outline orange}
}
}
# To plot all the cells, we go through the cell list and call plot for
# each cell.
proc plot_all {} {
global p
for {set j 0} {$j < $p(array_height)} {incr j} {
for {set i 0} {$i < $p(array_width)} {incr i} {
plot $j $i
}
}
for {set i 0} {$i < $p(array_width)} {incr i} {
plot_block $i
}
}
# We mark a gas cell by setting its type in its cell entry to "marked_n" where
# n is an integer that increments every time we call this routine. Thus all
# marked cells are marked with a unique type, which allows us to track them
# individually. This routine is designed to be called from a mouse click event at
# position (x,y) in canvas widget w. So we use (x,y) to find the cell rectangle
# the user clicked on. We set this cell's type to "marked_n" if it is "unmarked"
# or to "unmarked" if it is already marked. This is the routine that gets called
# when we click on a gas cell in the display.
proc mark_cell {x y w} {
if {![winfo exists $w]} {return}
global cells p
update
set name [$w find closest $x $y]
set index 0
for {set j 0} {$j < $p(array_height)} {incr j} {
for {set i 0} {$i < $p(array_width)} {incr i} {
if {[get_name $j $i] == $name} {
switch -glob [get_type $j $i] {
"unmarked" {
set_type $j $i "mark_[incr p(mark_num)]"
track $j $i
}
"mark*" {set_type $j $i "unmarked"}
default {set_type $j $i "unmarked"}
}
plot $j $i
update
break
}
}
}
}
# This routine allow us to switch a surface block from sand to water or some
# other material.
proc configure_block {x y w} {
if {![winfo exists $w]} {return}
global blocks p
update
set name [$w find closest $x $y]
set index 0
for {set i 0} {$i < $p(array_width)} {incr i} {
if {[get_block_name $i] == $name} {
set type [get_block_type $i]
switch -glob $type {
"sand" {set_block_type $i "water"}
"water" {set_block_type $i "sand"}
}
plot_block $i
update
break
}
}
}
# We reset all the cells by setting their temperatures to the initial temperature
# and their types to "unmarked". We set all surface blocks to "sand".
proc reset_cells {} {
global cells p
for {set j 0} {$j < $p(array_height)} {incr j} {
for {set i 0} {$i < $p(array_width)} {incr i} {
set_temp $j $i $p(T_initial)
set_type $j $i "unmarked"
set_wvc $j $i "0.0"
set_cdc $j $i "0.0"
set_sfc $j $i "0.0"
set_rdc $j $i "0.0"
}
}
for {set i 0} {$i < $p(array_width)} {incr i} {
set_block_temp $i $p(T_initial)
set_block_type $i "sand"
}
set p(counter) 0
plot_all
update
}
# We display a window with all the program parameters so the user can change them.
proc configure {} {
global p
set w .configure
if {[winfo exists $w]} {destroy $w}
toplevel $w
wm title $w "Circulating Cells, Configuration Array"
frame $w.f
pack $w.f -side left -fill x
foreach a $p(config_names) {
set b [string tolower $a]
label $w.f.l$b -text $a -anchor w -width 20
entry $w.f.e$b -textvariable p($a) -relief sunken -bd 1 -width 30
grid $w.f.l$b $w.f.e$b -sticky news
}
foreach a {cell_mass_kg cell_height_m p_top_kPa p_bottom_kPa} {
set b [string tolower $a]
label $w.f.l$b -text $a -anchor w -width 20 -fg blue
label $w.f.e$b -textvariable p($a) -width 30 -fg blue
grid $w.f.l$b $w.f.e$b -sticky news
}
}
# We open the data window, or bring it to the foreground, when we press
# the data button.
proc data {} {
global p
set width 60
set height 30
set w $p(data_window)
if {[winfo exists $w]} {
raise $w
return $w.text
}
toplevel $w
wm title $w "$p(name), Data Window"
frame $w.f
pack $w.f -side top -fill x
button $w.f.report -command report -text "Report"
pack $w.f.report
label $w.f.description -textvariable p(data_description) -fg blue
pack $w.f.description
frame $w.t
pack $w.t -side top -expand yes -fill both
set t [text $w.t.text -relief sunken -bd 1 \
-yscrollcommand "$w.t.vsb set" \
-xscrollcommand "$w.t.hsb set" \
-setgrid 1 -height $height -width $width -border 2 \
-wrap none -tabs "0.25i left" -undo 1]
scrollbar $w.t.vsb -orient vertical -command "$t yview"
scrollbar $w.t.hsb -orient horizontal -command "$t xview"
pack $w.t.vsb -side right -fill y
pack $w.t.hsb -side bottom -fill x
pack $t -expand yes -fill both
set p(text) $t
report
return $t
}
# This procedure saves the cell array to a file. It writes in comments
# the originator of the file and a few other details. The comment lines
# begin with a hash character.
proc save {{fn ""}} {
global cells p
# Stop the simulation in preparation for the save.
set p(control) "Stop"
# Allow the user to specify an array file.
if {$fn == ""} {
set fn [tk_getSaveFile -initialfile "cells.txt"]
}
if {$fn == ""} {return ""}
set f [open $fn w]
# Write some metadata to the file in comment lines.
puts $f "# Cell Array Data"
puts $f "# Created By: $p(name)"
foreach a $p(config_names) {
puts $f "# $a = \"$p($a)\""
}
# Write surface blocks to the file, each on a separate line.
for {set i 0} {$i < $p(array_width)} {incr i} {
set cell_line "block $i "
foreach a {type temp} {
append cell_line "[get_block_$a $i] "
}
puts $f $cell_line
}
# Write gas cells to the file, each on a separate line.
for {set j 0} {$j < $p(array_height)} {incr j} {
for {set i 0} {$i < $p(array_width)} {incr i} {
set cell_line "$j $i "
foreach a {type temp wvc cdc sfc rdc} {
append cell_line "[get_$a $j $i] "
}
puts $f $cell_line
}
}
# Close file and returen.
close $f
return $fn
}
# This procedure loads the cell array from a file. It also sets the
# cell counter to zero, on the assumption that by loading a new array
# you are starting afresh. The cell array file can contain comment lines
# that start with a hash character. If we pass a file name to the procedure,
# it uses that file. Otherwise it opens a browser to ask for a file.
proc load {{fn ""}} {
global cells p
# Stop the simulation in preparation for the load.
set p(control) "Stop"
# Allow the user to select an array file.
if {$fn == ""} {set fn [tk_getOpenFile]}
if {![file exists $fn]} {return ""}
# Read file contents.
set new_cells "\n"
set f [open $fn r]
append new_cells [read $f]
close $f
# From the comments, extract the version number
if {![regexp {Version ([^ ]+)\n} $new_cells match version]} {
set version 0.0
}
# Remove comment lines
regsub -all {\n[ \t]*#[^\n]*} $new_cells "" new_cells
# Split the cell array string into individual cells
# and over-write the old simulation array.
set new_cells [split [string trim $new_cells] \n]
# Transfer cells to our existing array list. Our interpretation
# of the cell strings depends upon the version of the program
# that wrote the settings file.
foreach c $new_cells {
set j [lindex $c 0]
set i [lindex $c 1]
if {$j == "block"} {
set_block_type $i [lindex $c 2]
set_block_temp $i [lindex $c 3]
} {
if {$version >= 7.0} {
set_type $j $i [lindex $c 2]
set_temp $j $i [lindex $c 3]
set_wvc $j $i [lindex $c 4]
set_cdc $j $i [lindex $c 5]
if {$version >= 10.0} {
set_sfc $j $i [lindex $c 6]
set_rdc $j $i [lindex $c 7]
} {
set_sfc $j $i "0.0"
set_rdc $j $i "0.0"
}
} {
set_temp $j $i [lindex $c 2]
set_type $j $i [lindex $c 4]
set_wvc $j $i "0.0"
set_cdc $j $i "0.0"
set_sfc $j $i "0.0"
set_rdc $j $i "0.0"
}
}
}
# Set the cell-counter to zero and exit.
set p(counter) 0
return $fn
}
# circulate performs rotation of a block of four cells by a quarter-turn if the
# the rotation will do enough work to accelerate the cells to a speed sufficient
# to bring about the rotation in the expected time that elapses between attempts
# to rotate these cells. We relate program iterations to time in our Simulation
# Time post. In Impetus for Circulation, we calculate the time available and
# relate it to the required energy. We call the work done by the rotation the
# impetus for circulation, and we explore in detail the various ways we can calculate
# the impetus in posts like Buoyancy Work, Expansion Work, and Impetus by Enthalpy.
# In the end, we find that the drop in the average temperature of the cells,
# multiplied by Cp, gives us the work available per kilogram of gas in the four
# cells. At the top of this program, we calculate the impetus threshold, which is
# the minimum work per kilogram of gas that is adequate to bring about the
# rotation.
proc circulate {j i} {
global p cells
# Then number of cells in the block.
set num_cells 4
# We will use j0, j1, i0, i1 for the positions within the block.
set j0 $j
set j1 [expr $j + 1]
set i0 $i
set i1 [expr $i + 1]
# The four cells in the block are Tji where j is 0 for the bottom row, 1 for top
# and i is 0 for left column and 1 for right. Our first step is to determine the
# temperatures of the four cells.
set T00 [get_temp $j0 $i0]
set T10 [get_temp $j1 $i0]
set T11 [get_temp $j1 $i1]
set T01 [get_temp $j0 $i1]
# We calculate the factor by which temperature changes when a cell rises from row
# j to row j+1. We start with p_factor, which is the proportional change in pressure
# from j to j+1.
set p_factor [expr (1.0-$p(p_step_fraction)*($j+1))/(1.0-$p(p_step_fraction)*$j)]
set T_factor [expr pow($p_factor,1.0*$p(R_gas)/$p(Cp_gas))]
# We calculate the impetus for circulation, which is the drop in the average
# temperature of the four cells, multiplied by Cp. The units are Joule per
# kilogram. There are two values of impetus: one for clockwise and one for
# anti-clockwise rotation. This calculation is within a few percent of the one
# we describe in our Buoyancy Work and Expansion Work posts, followed by
# Impetus Dissected. To connect the two calculations, consider the difference
# between the initial temperature of the rising cell and the final temperature
# of the falling cell. This is the delta-T in the Buoyancy Work diagram. If
# we calculate this here, multiply by mg (which is p_step), divide by p (which
# is p_bottom - j*p_step), multiply by R_gas, and divide by num_cells, we get
# an answer within a few percent of the one we obtain below.
set cw_impetus [expr $p(Cp_gas) \
* ($T00*(1.0-$T_factor) + $T11*(1.0-1.0/$T_factor)) \
/ $num_cells]
set acw_impetus [expr $p(Cp_gas) \
* ($T01*(1.0-$T_factor) + $T10*(1.0-1.0/$T_factor)) \
/ $num_cells]
# We rotate in the most favorable direction, or not at all. When we rotate, we
# multiply the temperature of rising cells by T_factor and divide the temperature
# of falling cells by the same factor. If the circulation takes place, we add the
# impetus for circulation back into the gas in the form of visous heat. In Work
# by Circulation, we showed that failing to account for this viscous dissipation
# resulted in 15% of the heat we put into our planetary greenhouse disappearing.
set rotation 0
if {($cw_impetus >= $acw_impetus) && ($cw_impetus >= $p(impetus_threshold))} {
set viscous_heating [expr $cw_impetus / $p(Cp_gas)]
set_temp $j0 $i0 [expr $T00 * $T_factor + $viscous_heating]
set_temp $j0 $i1 [expr $T01 + $viscous_heating]
set_temp $j1 $i1 [expr $T11 / $T_factor + $viscous_heating]
set_temp $j1 $i0 [expr $T10 + $viscous_heating]
swap $j0 $i0 $j0 $i1
swap $j0 $i1 $j1 $i1
swap $j1 $i1 $j1 $i0
set rotation 1
}
if {($acw_impetus > $cw_impetus) && ($acw_impetus >= $p(impetus_threshold))} {
set viscous_heating [expr $acw_impetus / $p(Cp_gas)]
set_temp $j0 $i0 [expr $T00 + $viscous_heating]
set_temp $j0 $i1 [expr $T01 * $T_factor + $viscous_heating]
set_temp $j1 $i1 [expr $T11 + $viscous_heating]
set_temp $j1 $i0 [expr $T10 / $T_factor + $viscous_heating]
swap $j0 $i0 $j1 $i0
swap $j1 $i0 $j1 $i1
swap $j1 $i1 $j0 $i1
set rotation -1
}
# We indicate which direction we rotated.
return $rotation
}
# The penetration procedure determines the fraction of the sun's heat that
# penetrates through the gas cells to the surface blocks. We obtain cloud
# depth in millimeters from the total mass of cloud water droplets in grams
# per square meter, then apply our exponential penetration depth. We assume
# that snow and rain are too concentrated to reflect a significant fraction of
# sunlight. It is the refraction by a multitude of droplets that reflects
# sunlight, not absorption by water.
proc penetration {i} {
global p
set sum 0.0
for {set j 0} {$j < $p(array_height)} {incr j} {
set sum [expr $sum + [get_cdc $j $i]]
}
set cloud_depth [expr $sum*$p(cell_mass_kg)/1000.0]
return [expr exp(-$cloud_depth/$p(Lc_water))]
}
# The heat routine controls how we introduce heat into the array and
# how we take it out. We call the routine only every heating_interval
# iterations so as to save execution time.
proc heat {} {
global p
set j_max [expr $p(array_height) - 1]
switch $p(time_control) {
"Day" {
set p(solar_time) $p(mid_day_hr)
set p(Q_current) $p(Q_sun)
}
"Night" {
set p(solar_time) $p(mid_night_hr)
set p(Q_current) 0.0
}
"Cycle" {
set p(solar_time) [format %.2f [expr \
fmod(1.0 * $p(counter) / $p(iterations_per_s) / $p(s_per_hr), \
$p(day_length_hr))]]
if {[expr abs(2.0*($p(solar_time)-$p(mid_day_hr))/$p(day_length_hr))] \
< $p(daylight_fraction)} {
set p(Q_current) $p(Q_sun)
} {
set p(Q_current) 0.0
}
}
}
# We are going to make a note of the total heat radiated into space
# by the atmosphere, so here we set a variable to zero that will
# add it all up, and another for the long-wave radiation descending
# from the atmospheric clouds and gas.
set Q_upwelling_sum 0.0
set Q_downwelling_sum 0.0
for {set i 0} {$i < $p(array_width)} {incr i} {
# We are going to calculate the heat gained by the
# surface block and each cell above it. Negative values
# represent heat lost. We use the variable Q_n to store
# the power per square meter gained by cell n, and Q_b
# for the surface block. The first thing we do is set
# these to zero. Later, we multiply the net power per
# square meter by the time between calls to this heating
# routine, and so obtain the total heat gained.
set Q_b 0.0
for {set j 0} {$j < $p(array_height)} {incr j} {
set Q_$j 0.0
}
# Determine the time for which we will calculate
# heat transfer.
set t [expr $p(heating_interval) / $p(iterations_per_s)]
# The heat arriving at the surface block from the sun.
# For this we call the penetration routine, which calculates the
# fraction of sunlight that makes it through the clouds.
set Q_b [expr $Q_b + $p(Q_current) * [penetration $i]]
# Obtain the block temperature and heat capacity
# per square meter.
set Tb [get_block_temp $i]
set Cb [set p(heat_capacity_[get_block_type $i])]
# Obtain the surface gas cell temperature and capacity per
# square meter.
set Tc [get_temp 0 $i]
set Cc [expr $p(cell_mass_kg) * $p(Cp_gas)]
# The convection transfer rate to the gas cell above.
if {$Tb > $Tc} {
set Q [expr ($Tb - $Tc) * $p(convection_coeff)]
set Q_0 [expr $Q_0 + $Q]
set Q_b [expr $Q_b - $Q]
}
# The evaporation of water per square meter of the surface and the
# accompanying loss of latent heat. We make sure that the evaporation
# does not cause the concentration of water vapor in the cell to
# exceed saturation. We determine the heat carried away from the
# latent heat of evaporation in W/m2.
if {[get_block_type $i] == "water"} {
set Xs [saturation_concentration $Tc 0]
set Xc [get_wvc 0 $i]
if {$Xc < $Xs} {
set W [evaporation_rate $Tb [get_wvc 0 $i]]
if {$Xc + $W*$t/$p(cell_mass_kg) > $Xs} {
set W [expr ($Xs-$Xc)*$p(cell_mass_kg)/$t]
}
set Q_b [expr $Q_b - $W * $p(Qe_water)]
set_wvc 0 $i [expr [get_wvc 0 $i] + $W * $t / $p(cell_mass_kg)]
}
}
# We now ascend from the surface block up to the tropopause
# calculating the upward long-wave radiation power, on the
# assumption that it goes verticaly. We start with the heat
# radiated by the surface block, which is always a black body.
# Because the dry gas is opaque to some wavelengths and
# transparent to others, we separate the heat into that which
# is in the transparent wavelengths and that which is in the
# opaque wavelengths.
set Q [expr pow($Tb,4.0) * $p(stefan_const)]
set Q_b [expr $Q_b - $Q]
set Q_t [expr $Q * $p(transparency_fraction)]
set Q_o [expr $Q - $Q_t]
# With these initial value of Q_t and Q_o, we step through the
# cells above. If the cell is a cloud, it absorbs all upward
# radiation and radiates upwards as a black body. Otherwise it
# absorbs only the opaque wavelength radiation.
for {set j 0} {$j < $p(array_height)} {incr j} {
set T [get_temp $j $i]
if {[get_cdc $j $i]+[get_sfc $j $i]+[get_rdc $j $i] >= $p(wc_opaque)} {
set Q_$j [expr [set Q_$j] + $Q_t + $Q_o]
set Q [expr pow($T,4.0) * $p(stefan_const)]
set Q_t [expr $Q * $p(transparency_fraction)]
set Q_o [expr $Q - $Q_t]
set Q_$j [expr [set Q_$j] - $Q]
} {
set Q_$j [expr [set Q_$j] + $Q_o]
set Q_o [expr pow($T,4.0) * $p(stefan_const) \
* (1-$p(transparency_fraction))]
set Q_$j [expr [set Q_$j] - $Q_o]
}
}
# At this point Q_o + Q_t are the total heat radiated by the surface
# and clouds upwards.
set Q_upwelling_sum [expr $Q_upwelling_sum + $Q_o + $Q_t]
# Now we start at the top row, the tropopause, and work our way down,
# calculating the downward long-wave power, using Q_o and Q_t in the
# same way, but in the opposite direction. We set both to zero to start
# with because no long-wave radiation arrives from outer space.
set Q_o 0.0
set Q_t 0.0
for {set j [expr $p(array_height)-1]} {$j >= 0} {incr j -1} {
set T [get_temp $j $i]
if {[get_cdc $j $i]+[get_sfc $j $i]+[get_rdc $j $i] >= $p(wc_opaque)} {
set Q_$j [expr [set Q_$j] + $Q_t + $Q_o]
set Q [expr pow($T,4.0) * $p(stefan_const)]
set Q_t [expr $Q * $p(transparency_fraction)]
set Q_o [expr $Q - $Q_t]
set Q_$j [expr [set Q_$j] - $Q]
} {
set Q_$j [expr [set Q_$j] + $Q_o]
set Q_o [expr pow($T,4.0) * $p(stefan_const) \
* (1-$p(transparency_fraction))]
set Q_$j [expr [set Q_$j] - $Q_o]
}
}
# All the downward going radiation at the bottom of the atmosphere is
# absorbed by the surface block.
set Q_b [expr $Q_b + $Q_o + $Q_t]
# At this point, Q_o + Q_t are the total heat radiated to the surface
# by the clouds and atmosphere so we add them to the downwelling sum.
set Q_downwelling_sum [expr $Q_downwelling_sum + $Q_o + $Q_t]
# Adjust the surface block temperature.
set_block_temp $i [expr [get_block_temp $i] + $Q_b * $t / $Cb]
# Adjust the cell temperatures.
for {set j 0} {$j < $p(array_height)} {incr j} {
set_temp $j $i [expr [get_temp $j $i] + [set Q_$j] * $t / $Cc]
}
if {!$p(accelerate)} {
plot 0 $i
plot_block $i
}
}
# We record the average radiation into space and the average radiation
# received from the atmosphere.
set p(Q_downwelling) [expr $Q_downwelling_sum / $p(array_width)]
set p(Q_upwelling) [expr $Q_upwelling_sum / $p(array_width)]
return 1
}
# check_cells will stop the simulation by setting control to "Stop"
# when certain user-specified termination criteria are met. This is the
# place to put termination code.
proc check {} {
global p
set j_max [expr $p(array_height) - 1]
set i_max [expr $p(array_width) - 1]
# begin code-pasting area.
# end code-pasting area.
}
# report writes array properties to the data window if it exists.
proc report {} {
global p
if {![winfo exists $p(data_window)]} {return 0}
# Set up the new line and the description string.
set line ""
set description ""
# The time in hours.
set sim_time [format %.2f [expr 1.0*$p(counter)/$p(iterations_per_s)/$p(s_per_hr)]]
append line "$sim_time "
append description "time (hr) "
# The average sand block temperature.
set sum 0.0
set count 0
for {set i 0} {$i < $p(array_width)} {incr i} {
if {[get_block_type $i] != "water"} {
set sum [expr $sum + [get_block_temp $i]]
incr count
}
}
if {$count > 0} {
append line "[format %.2f [expr $sum/$count]] "
} {
append line "-1 "
}
append description "sand (K) "
# The average water block temperature.
set sum 0.0
set count 0
for {set i 0} {$i < $p(array_width)} {incr i} {
if {[get_block_type $i] == "water"} {
set sum [expr $sum + [get_block_temp $i]]
incr count
}
}
if {$count > 0} {
append line "[format %.2f [expr $sum/$count]] "
} {
append line "-1 "
}
append description "water (K) "
# Average surface gas cell.
append line "[get_row_temp 0] "
append description "surface (K) "
# Average tropopause temperature.
append line "[get_row_temp [expr $p(array_height)-1]] "
append description "tropopause (K) "
# Average cloud depth in millimeters, including snow and droplets,
# but not rain or water vapor.
set sum 0.0
for {set i 0} {$i < $p(array_width)} {incr i} {
for {set j 0} {$j < $p(array_height)} {incr j} {
set sum [expr $sum + [get_cdc $j $i]]
}
}
set cloud_depth [expr $sum*$p(cell_mass_kg)/1000.0/$p(array_width)]
append line "[format %.2f $cloud_depth] "
append description "clouds (mm) "
# The incoming solar heat, before reflection by clouds.
append line "[format %.2f $p(Q_current)] "
append description "solar (W/m2) "
# Average heat from the sun in Watt per square meter.
set sum 0.0
for {set i 0} {$i < $p(array_width)} {incr i} {
set sum [expr $sum + [penetration $i] * $p(Q_current)]
}
append line "[format %.2f [expr $sum/$p(array_width)]] "
append description "penetrating (W/m2) "
# Average upwelling radiation.
append line "[format %.2f $p(Q_upwelling)] "
append description "up (W/m2) "
# Average downwelling radiation.
append line "[format %.2f $p(Q_downwelling)] "
append description "down (W/m2) "
# We write our total line to the text window and put the description
# in the data window label.
if {1} {
set p(data_description) $description
$p(text) insert end "$line\n"
$p(text) see end
}
# Here are some other reports we use, with if statements to turn them
# off. One gets the row temperatures and the other the row altitudes.
# Set the 0 to 1 to turn on the code, and 1 to 0 to turn it off.
if {0} {
for {set j 0} {$j < $p(array_height)} {incr j} {
$p(text) insert end "[get_row_altitude $j] "
}
$p(text) insert end "\n"
}
if {0} {
for {set j 0} {$j < $p(array_height)} {incr j} {
$p(text) insert end "[get_row_temp $j] "
}
$p(text) insert end "\n"
}
# The following lines we use to increase the solar power and record
# the state of the array at fixed intervals of one hundred times
# the reporting interval.
if {0} {
if {($p(counter) > 0) \
&& (fmod($p(counter),100*$p(report_interval)) == 0.0)} {
save "Cells_$p(Q_sun)\.txt"
set p(control) "Run"
set f [open "Data.txt" a]
puts $f $line
close $f
set p(Q_sun) [expr $p(Q_sun) + 10.0]
}
}
return $line
}
# track writes the properties of a cell to the text window, should
# it be open.
proc track {j i} {
global cells p
if {!$p(tracking)} {return 0}
if {[winfo exists $p(data_window)]} {
$p(text) insert end "[get_type $j $i] $j $i "
set Tc [get_temp $j $i]
set Xc [get_wvc $j $i]
set Dc [get_cdc $j $i]
set Fc [get_sfc $j $i]
set Rc [get_rdc $j $i]
set Xs [saturation_concentration $Tc $j]
set Ts [dew_point $Tc $Xs $j]
set pc [expr $p(p_bottom_kPa) - $j * $p(p_step_kPa)]
foreach a {Tc Xc Dc Fc Rc Xs Ts pc} {
$p(text) insert end "[format %.2f [set $a]] "
}
$p(text) insert end "\n"
$p(text) see end
}
}
# We start an infinite loop whose actions are controlled by the buttons
# at the top of the display. The loop pays careful attention to how often
# it responds to user input and how often it re-plots the cells. It runs
# in a slow or a fast mode, depending upon the state of the accelerate
# flag. Each run through the loop calls the array heating routine, selects
# a random block of four cells, and rotates it if it's ready to rotate, and
# checks the water balance within the cell.
while {1} {
if {[catch {winfo exists .}]} {exit}
if {$p(control) == "Reset"} {
reset_cells
set p(control) "Stop"
}
if {$p(control) == "Stop"} {
plot_all
update
after 100
continue
}
# We perform various operations with a frequency set by their "interval".
if {[expr fmod($p(counter),$p(plot_interval))] == 0} {plot_all}
if {([expr fmod($p(counter),$p(report_interval))] == 0) \
&& $p(reporting)} {report}
if {[expr fmod($p(counter),$p(check_interval))] == 0} {check}
# We update the display on every iteration unless we are operating in
# accelerated mode, in which case we update less often.
if {[expr fmod($p(counter),$p(update_interval))] == 0} {
update
} else {
if {!$p(accelerate)} {
update
}
}
# We heat the cells every heating_interval to save execution time.
if {[expr fmod($p(counter),$p(heating_interval))] == 0} {heat}
# We pick a block of four cells at random, specifying it with its bottom-left
# cell.
set j [expr round(rand()*($p(array_height)-1)-0.5)]
set i [expr round(rand()*($p(array_width)-1)-0.5)]
# We apply circulation.
set rotation [circulate $j $i]
# We apply water balance to each cell.
foreach jj "$j [expr $j+1]" {
foreach ii "$i [expr $i+1]" {
water_balance $jj $ii
}
}
# We plot the cells in the block, and those beneath them, if they
# exist, because they may have received rain or snow.
foreach jj "[expr $j - 1] $j [expr $j+1]" {
if {$jj >= 0} {
foreach ii "$i [expr $i+1]" {
plot $jj $ii
}
}
}
# We provide tracking for marked gas cells each time we check them.
if {$p(tracking)} {
foreach jj "$j [expr $j+1]" {
foreach ii "$i [expr $i+1]" {
if {[get_type $jj $ii] != "unmarked"} {track $jj $ii}
}
}
}
# We increment the iteration counter.
incr p(counter)
}