Menu
▾
▴
perl source code
!/usr/bin/perl -w
use strict;
use Statistics::Descriptive;
File: run_container.pl
Description: Runs container modeling program for specified experiments. Constructs parameter file based on a template, then runs program. Runs plots.
Input: Container program, specify parameters here, input met data
Output: Container program output files
Dan Steinhoff
13 March 2013
USER SPECIFICATIONS
Location
my @cities = qw/ Veracruz RioBlanco Puebla /; # List of cities
my @sublocs = qw/ Bachilleres RioBlanco Carmen /; # List of sublocations within cities
my @lats = (19.1983, 18.8367, 19.0375); # Site latitudes
my @lons = (-96.1350, -97.1511, -98.2397); # Site longitudes
my @elevs = (56.08, 1270.2, 2114.70); # Site elevations
Containers
my @cshapes = qw/ ROUND ROUND ROUND /; # Shape of containers ('ROUND', 'RECT', 'TIRE')
my @radii1_t = (0.082, 0.131, 0.298); # Top container lengths 1 (in m)
my @radii2_t = (0.082, 0.131, 0.298); # Top container lengths 2 (in m)
my @radii1_b = (0.082, 0.131, 0.298); # Body container lengths 1 (in m)
my @radii2_b = (0.082, 0.131, 0.298); # Body container lengths 2 (in m)
my @heights = (0.191, 0.368, 0.921); # Container heights (in m)
my @conducs = (0.30, 0.50, 0.50); # Thermal conductivity of containers (in W m-1 K-1)
my @thicks = (0.0018, 0.0023, 0.0048); # Thickness of containers (in m)
my @names_cont = qw/ Coffee Bucket Barrel /; # Names of containers
Reflectivities (colors)
my @refls = (0.1, 0.5, 0.9); # Reflectivity coefficients of containers - depends on number of containers
my @names_refl = qw/ Black Gray White /; # Names of reflectivities (colors) - depends on number of containers
Shade
my @shades = (0.00, 0.50, 1.00); # Fraction of shade
my @names_shade = qw/ NoShade HalfShade FullShade /; # Names of shade conditions
Clouds (if not provided)
my @clouds = ([0.0,0.0,0.0,0.0,0.0,0.0], [0.33,0.33,0.33,0.1,0.1,0.1], [0.66,0.66,0.66,0.33,0.33,0.33]); # Inner index is low, middle, and high clouds, for day and night.
my @names_cloud = qw/ Clear PartlyCloudy Overcast /; # Names of cloud cover conditions (if specified)
Water, ground, and cloud conditions
my @whs_mul = (0.0, 0.0, 0.0); # Multiplier of initial water height / container height ratio - depends on number of containers
my @tws = (25.00, 20.00, 18.00); # Initial water temperatures (in deg C) - depends on number of locations
my @tgs = (26.86, 18.82, 18.16); # Initial ground temperatures - depends on number of locations
Flags
my $cflag = "F"; # Whether cloud fraction data is input (T/F) | leer archivo externo para utilizar datos de nubes Si (T)/no(F)
my $sflag = "F"; # Whether soil temperature data is input (T/F) | leer archivo externo para para utilizar datos de temperatura del suelo Si (T)/no(F)
my $mfill = "T"; # Manual fill flag (T/F) | llenar contenedores manualmente Si (T)/no(F)
my $fill_intv = 12; # Manual fill interval, in hours
Misc.
my $time_step = 60.00; # Model time step in seconds
my $input_intv = 3600.00; # Time interval of input data in seconds
my $miss = -9999.00; # Missing value of input data
my $plot = 1; # Whether to run plotting routines or not (1/0)
my $path = "/home/roger/Dropbox/watchem/"; # Path to container.exe file | Ruta al archivo ejecutable container.exe
END USER SPECIFICATIONS
check line 87 to declare meteo data source files.
Welcome message]
print "##############################\n";
print "#### WELCOME TO WHATCH'EM ####\n";
print "##############################\n";
print "# Starting...\n";
Get array sizes for use throughout program
my $exp_cont_size = scalar (@radii1_t);
my $exp_shade_size = scalar (@shades);
my $exp_refl_size = scalar (@refls);
my $exp_cloud_size = scalar (@clouds);
experimento 1
Open master statistics files for writing | ruta a archivos maestros de estadisticos climaticos para escribir
open MTOUTFILE, ">/home/roger/Dropbox/watchem/stats/stats_temperature.txt" or die "Cannot open master temperature stats file to write: $!";
open MWOUTFILE, ">/home/roger/Dropbox/watchem/stats/stats_waterheight.txt" or die "Cannot open master water height stats file to write: $!";
if ($cflag eq 'F')
{
printf MTOUTFILE "%-30s %-10s %-15s %-10s %-15s %-6s %-6s %-6s %-6s %-6s %-6s %-6s %-6s\n", "CITY-HOBO", "CONTAINER", "REFLS", "SHADE", "CLOUDS", "TMEAN", "TMAX", "TMIN", "CMEAN", "CMAX", "CMIN", "TAMEAN", "EVAP";
printf MWOUTFILE "%-30s %-10s %-15s %-10s %-15s %-15s\n", "CITY-HOBO", "CONTAINER", "REFLS", "SHADE", "CLOUDS", "WATER HEIGHT";
}
else
{
printf MTOUTFILE "%-30s %-10s %-15s %-10s %-6s %-6s %-6s %-6s %-6s %-6s %-6s %-6s\n", "CITY-HOBO", "CONTAINER", "REFLS", "SHADE", "TMEAN", "TMAX", "TMIN", "CMEAN", "CMAX", "CMIN", "TAMEAN", "EVAP";
printf MWOUTFILE "%-30s %-10s %-15s %-10s %-15s\n", "CITY-HOBO", "CONTAINER", "REFLS", "SHADE", "WATER HEIGHT";
}
Loop through each city
my @sublocs_t = @sublocs;
foreach my $city (@cities)
{
print "# $city\n";
# Get the city-specific variables
my $subloc = shift (@sublocs_t);
my $input_file = "/home/roger/Dropbox/watchem/input/container_input_$city-$subloc.txt";
my $output_file_base = "container_output_$city-$subloc";
if (!(-d "output/$city-$subloc")) { system "mkdir output/$city-$subloc" }
my $lat = sprintf "%.2f", shift (@lats);
my $lon = sprintf "%.2f", shift (@lons);
my $elev = sprintf "%.2f", shift (@elevs);
my $tw = sprintf "%.2f", shift (@tws);
my $tg = sprintf "%.2f", shift (@tgs);
my $t_diff = $input_intv/$time_step;
$time_step = sprintf "%.2f", $time_step;
$t_diff = sprintf "%.2f", $t_diff;
# Loop through each container for (my $i = 0; $i < $exp_cont_size; $i++) { my $cshape = $cshapes[$i]; my $radius1_t = sprintf "%.4f", $radii1_t[$i]; my $radius2_t = sprintf "%.4f", $radii2_t[$i]; my $radius1_b = sprintf "%.4f", $radii1_b[$i]; my $radius2_b = sprintf "%.4f", $radii2_b[$i]; my $height = sprintf "%.4f", $heights[$i]; my $wh = sprintf "%.4f", $height*$whs_mul[$i]; my $name_cont = $names_cont[$i]; print "## $name_cont\n"; # Loop through each reflectivity for (my $j = 0; $j < $exp_refl_size; $j++) { my $name_refl = $names_refl[$j]; print "### $name_refl\n"; my $refl = sprintf "%.2f", $refls[$j]; my $conduc = sprintf "%.2f", $conducs[$j]; my $thick = sprintf "%.4f", $thicks[$j]; # Loop through each shade for (my $k = 0; $k < $exp_shade_size; $k++) { my $name_shade = $names_shade[$k]; print "#### $name_shade\n"; my $shade = sprintf "%.2f", $shades[$k]; # Loop through each cloud (if not provided) if ($cflag eq 'F') { for (my $l = 0; $l < $exp_cloud_size; $l++) { my $name_cloud = $names_cloud[$l]; print "##### $name_cloud\n"; my $output_file = $output_file_base."-$name_cont-$name_refl-$name_shade-$name_cloud.txt"; my @clouds_exp = @{$clouds[$l]}; my $cld = sprintf "%.2f", $clouds_exp[0]; my $cmd = sprintf "%.2f", $clouds_exp[1]; my $chd = sprintf "%.2f", $clouds_exp[2]; my $cln = sprintf "%.2f", $clouds_exp[3]; my $cmn = sprintf "%.2f", $clouds_exp[4]; my $chn = sprintf "%.2f", $clouds_exp[5]; # Construct parameter file print "##### Constructing parameter file...\n"; open OUTFILE, ">param_file.txt"; print OUTFILE "$input_file\n"; print OUTFILE "$city\n"; print OUTFILE "$subloc\n"; print OUTFILE "$lat\n"; print OUTFILE "$lon\n"; print OUTFILE "$elev\n"; print OUTFILE "$cshape\n"; print OUTFILE "$radius1_t\n"; print OUTFILE "$radius2_t\n"; print OUTFILE "$radius1_b\n"; print OUTFILE "$radius2_b\n"; print OUTFILE "$height\n"; print OUTFILE "$refl\n"; print OUTFILE "$conduc\n"; print OUTFILE "$thick\n"; print OUTFILE "$wh\n"; print OUTFILE "$tw\n"; print OUTFILE "$tg\n"; print OUTFILE "$time_step\n"; print OUTFILE "$t_diff\n"; print OUTFILE "$shade\n"; print OUTFILE "$cflag\n"; print OUTFILE "$sflag\n"; print OUTFILE "$mfill\n"; print OUTFILE "$fill_intv\n"; print OUTFILE "$output_file\n"; print OUTFILE "$cld\n"; print OUTFILE "$cmd\n"; print OUTFILE "$chd\n"; print OUTFILE "$cln\n"; print OUTFILE "$cmn\n"; print OUTFILE "$chn\n"; close OUTFILE; # Run container modeling program print "##### Running container model..."; system "$path/container.exe";
print "How old are you?";
$age = <>;
print "WOW! You are $age years old!";
# Move output file to correct directory print "##### Writing output files...\n"; system "mv $output_file output/$city-$subloc"; # Run daily energy balance and basic stats my $einfile = "output/$city-$subloc/$output_file"; my $eoutfile = "output/$city-$subloc/energybal_avg_$city-$subloc-$name_cont-$name_refl-$name_shade-$name_cloud.txt"; my $soutfile = "output/$city-$subloc/stats_$city-$subloc-$name_cont-$name_refl-$name_shade-$name_cloud.txt"; &daily_energy_balance_stats($einfile, $eoutfile, $soutfile, "$city-$subloc", $name_cont, $name_refl, $name_shade, $cflag, $name_cloud); # Run first set of plots if ($plot) { print "##### Running temperature and energy balance plots...\n"; # Plotting - Temperature variables (tw, ta, tg) system "ncl \'infile=\"output/$city-$subloc/$output_file\"\' \'outfile=\"plots/output/temperature/$city-${subloc}-$name_cont-$name_refl-$name_shade-$name_cloud\"\' \'container=\"$name_cont\"\' \'refl=\"$name_refl\"\' \'shade=\"$name_shade\"\' \'clouds=\"$name_cloud\"\' \'site=\"$city-$subloc\"\' plots/plot_t.ncl"; # Plotting - Average energy balance - Water system "ncl \'infile=\"output/$city-$subloc/energybal_avg_$city-$subloc-$name_cont-$name_refl-$name_shade-$name_cloud.txt\"\' \'outfile1=\"plots/output/energy_balance/$city-$subloc-$name_cont-$name_refl-$name_shade-water\"\' \'outfile2=\"plots/output/energy_balance/$city-$subloc-$name_cont-$name_refl-${name_shade}\"\' \'container=\"$name_cont\"\' \'refl=\"$name_refl\"\' \'shade=\"$name_shade\"\' \'clouds=\"$name_cloud\"\' \'site=\"$city-$subloc\"\' plots/plot_energy_balance_avg_water.ncl"; # Plotting - Average energy balance - Container system "ncl \'infile=\"output/$city-$subloc/energybal_avg_$city-$subloc-$name_cont-$name_refl-$name_shade-$name_cloud.txt\"\' \'outfile=\"plots/output/energy_balance/$city-$subloc-$name_cont-$name_refl-$name_shade-container\"\' \'container=\"$name_cont\"\' \'refl=\"$name_refl\"\' \'shade=\"$name_shade\"\' \'clouds=\"$name_cloud\"\' \'site=\"$city-$subloc\"\' plots/plot_energy_balance_avg_con.ncl"; # Plotting = Average daily water temperature system "ncl \'infile=\"output/$city-$subloc/energybal_avg_$city-$subloc-$name_cont-$name_refl-$name_shade-$name_cloud.txt\"\' \'outfile=\"plots/output/avg_temperature/$city-$subloc-$name_cont-$name_refl-$name_shade-avg\"\' \'container=\"$name_cont\"\' \'refl=\"$name_refl\"\' \'shade=\"$name_shade\"\' \'clouds=\"$name_cloud\"\' \'site=\"$city-$subloc\"\' plots/plot_t_avg.ncl"; } } # Plotting - Cloud experiments if ($plot) { print "##### Running cloud experiment plots...\n"; system "ncl \'infile1=\"output/$city-$subloc/container_output_$city-$subloc-$name_cont-$name_refl-$name_shade-Clear.txt\"\' \'infile2=\"output/$city-$subloc/container_output_$city-$subloc-$name_cont-$name_refl-$name_shade-PartlyCloudy.txt\"\' \'infile3=\"output/$city-$subloc/container_output_$city-$subloc-$name_cont-$name_refl-$name_shade-Overcast.txt\"\' \'outfile1=\"plots/output/comparison_clouds/temperature-$city-${subloc}-$name_cont-$name_refl-$name_shade-clouds\"\' \'outfile2=\"plots/output/comparison_clouds/waterheight-$city-${subloc}-$name_cont-$name_refl-$name_shade-clouds\"\' \'container=\"$name_cont\"\' \'refl=\"$name_refl\"\' \'shade=\"$name_shade\"\' \'site=\"$city-$subloc\"\' plots/plot_comparison_cloud.ncl"; } }
else
{
my $output_file = $output_file_base."-$name_cont-$name_refl-$name_shade.txt";
# Construct parameter file
print "#### Constructing parameter file...\n";
open OUTFILE, ">param_file.txt";
print OUTFILE "$input_file\n";
print OUTFILE "$city\n";
print OUTFILE "$subloc\n";
print OUTFILE "$lat\n";
print OUTFILE "$lon\n";
print OUTFILE "$elev\n";
print OUTFILE "$cshape\n";
print OUTFILE "$radius1_t\n";
print OUTFILE "$radius2_t\n";
print OUTFILE "$radius1_b\n";
print OUTFILE "$radius2_b\n";
print OUTFILE "$height\n";
print OUTFILE "$refl\n";
print OUTFILE "$conduc\n";
print OUTFILE "$thick\n";
print OUTFILE "$wh\n";
print OUTFILE "$tw\n";
print OUTFILE "$tg\n";
print OUTFILE "$time_step\n";
print OUTFILE "$t_diff\n";
print OUTFILE "$shade\n";
print OUTFILE "$cflag\n";
print OUTFILE "$sflag\n";
print OUTFILE "$mfill\n";
print OUTFILE "$fill_intv\n";
print OUTFILE "$output_file\n";
close OUTFILE;
# Run container modeling program print "#### Running container model..."; system "$path/container.exe"; # Move output file to correct directory print "#### Writing output files...\n"; system "mv $output_file output/$city-$subloc"; # Run daily energy balance and basic stats my $einfile = "output/$city-$subloc/$output_file"; my $eoutfile = "output/$city-$subloc/energybal_avg_$city-$subloc-$name_cont-$name_refl-$name_shade.txt"; my $soutfile = "output/$city-$subloc/stats_$city-$subloc-$name_cont-$name_refl-$name_shade.txt"; &daily_energy_balance_stats($einfile, $eoutfile, $soutfile, "$city-$subloc", $name_cont, $name_refl, $name_shade, $cflag); # Plotting routines if ($plot) { print "#### Running temperature and energy balance plots...\n"; # Plotting - Temperature variables (tw, ta, tg) system "ncl \'infile=\"output/$city-$subloc/$output_file\"\' \'outfile=\"plots/output/temperature/$city-${subloc}-$name_cont-$name_refl-$name_shade\"\' \'container=\"$name_cont\"\' \'refl=\"$name_refl\"\' \'shade=\"$name_shade\"\' \'clouds=\"Observed\"\' \'site=\"$city-$subloc\"\' plots/plot_t.ncl"; # Plotting - Average energy balance - Water system "ncl \'infile=\"output/$city-$subloc/energybal_avg_$city-$subloc-$name_cont-$name_refl-$name_shade.txt\"\' \'outfile1=\"plots/output/energy_balance/$city-$subloc-$name_cont-$name_refl-$name_shade-water\"\' \'outfile2=\"plots/output/energy_balance/$city-$subloc-$name_cont-$name_refl-$name_shade\"\' \'container=\"$name_cont\"\' \'refl=\"$name_refl\"\' \'shade=\"$name_shade\"\' \'clouds=\"Observed\"\' \'site=\"$city-$subloc\"\' plots/plot_energy_balance_avg_water.ncl"; # Plotting - Average energy balance - Container system "ncl \'infile=\"output/$city-$subloc/energybal_avg_$city-$subloc-$name_cont-$name_refl-$name_shade.txt\"\' \'outfile=\"plots/output/energy_balance/$city-$subloc-$name_cont-$name_refl-$name_shade-container\"\' \'container=\"$name_cont\"\' \'refl=\"$name_refl\"\' \'shade=\"$name_shade\"\' \'clouds=\"Observed\"\' \'site=\"$city-$subloc\"\' plots/plot_energy_balance_avg_con.ncl"; # Plotting = Average daily water temperature system "ncl \'infile=\"output/$city-$subloc/energybal_avg_$city-$subloc-$name_cont-$name_refl-$name_shade.txt\"\' \'outfile=\"plots/output/avg_temperature/$city-$subloc-$name_cont-$name_refl-$name_shade-avg\"\' \'container=\"$name_cont\"\' \'refl=\"$name_refl\"\' \'shade=\"$name_shade\"\' \'clouds=\"Observed\"\' \'site=\"$city-$subloc\"\' plots/plot_t_avg.ncl"; } } } } if ($plot) { # Plotting - Reflectivity experiments print "#### Running color experiment plots...\n"; for (my $k = 0; $k < $exp_shade_size; $k++) { my $name_shade = $names_shade[$k]; if ($cflag eq 'F') { for (my $l = 0; $l < $exp_cloud_size; $l++) { my $name_cloud = $names_cloud[$l]; system "ncl \'infile1=\"output/$city-$subloc/container_output_$city-$subloc-$name_cont-Black-$name_shade-$name_cloud.txt\"\' \'infile2=\"output/$city-$subloc/container_output_$city-$subloc-$name_cont-Gray-$name_shade-$name_cloud.txt\"\' \'infile3=\"output/$city-$subloc/container_output_$city-$subloc-$name_cont-White-$name_shade-$name_cloud.txt\"\' \'outfile1=\"plots/output/comparison_refls/temperature-$city-${subloc}-$name_cont-$name_shade-$name_cloud-refls\"\' \'outfile2=\"plots/output/comparison_refls/waterheight-$city-${subloc}-$name_cont-$name_shade-$name_cloud-refls\"\' \'container=\"$name_cont\"\' \'shade=\"$name_shade\"\' \'clouds=\"$name_cloud\"\' \'site=\"$city-$subloc\"\' plots/plot_comparison_refl.ncl"; } } else { system "ncl \'infile1=\"output/$city-$subloc/container_output_$city-$subloc-$name_cont-Black-$name_shade.txt\"\' \'infile2=\"output/$city-$subloc/container_output_$city-$subloc-$name_cont-Gray-$name_shade.txt\"\' \'infile3=\"output/$city-$subloc/container_output_$city-$subloc-$name_cont-White-$name_shade.txt\"\' \'outfile1=\"plots/output/comparison_refls/temperature-$city-${subloc}-$name_cont-$name_shade-refls\"\' \'outfile2=\"plots/output/comparison_refls/waterheight-$city-${subloc}-$name_cont-$name_shade-refls\"\' \'container=\"$name_cont\"\' \'shade=\"$name_shade\"\' \'clouds=\"Observed\"\' \'site=\"$city-$subloc\"\' plots/plot_comparison_refl.ncl"; } } # Plotting - Shade experiments print "### Running shade experiment plots...\n"; for (my $j = 0; $j < $exp_refl_size; $j++) { my $name_refl = $names_refl[$j]; if ($cflag eq 'F') { for (my $l = 0; $l < $exp_cloud_size; $l++) { my $name_cloud = $names_cloud[$l]; system "ncl \'infile1=\"output/$city-$subloc/container_output_$city-$subloc-$name_cont-$name_refl-NoShade-$name_cloud.txt\"\' \'infile2=\"output/$city-$subloc/container_output_$city-$subloc-$name_cont-$name_refl-HalfShade-$name_cloud.txt\"\' \'infile3=\"output/$city-$subloc/container_output_$city-$subloc-$name_cont-$name_refl-FullShade-$name_cloud.txt\"\' \'outfile1=\"plots/output/comparison_shade/temperature-$city-${subloc}-$name_cont-$name_refl-$name_cloud-shade\"\' \'outfile2=\"plots/output/comparison_shade/waterheight-$city-${subloc}-$name_cont-$name_refl-$name_cloud-shade\"\' \'container=\"$name_cont\"\' \'refl=\"$name_refl\"\' \'clouds=\"$name_cloud\"\' \'site=\"$city-$subloc\"\' plots/plot_comparison_shade.ncl"; } } else { system "ncl \'infile1=\"output/$city-$subloc/container_output_$city-$subloc-$name_cont-$name_refl-NoShade.txt\"\' \'infile2=\"output/$city-$subloc/container_output_$city-$subloc-$name_cont-$name_refl-HalfShade.txt\"\' \'infile3=\"output/$city-$subloc/container_output_$city-$subloc-$name_cont-$name_refl-FullShade.txt\"\' \'outfile1=\"plots/output/comparison_shade/temperature-$city-${subloc}-$name_cont-$name_refl-shade\"\' \'outfile2=\"plots/output/comparison_shade/waterheight-$city-${subloc}-$name_cont-$name_refl-shade\"\' \'container=\"$name_cont\"\' \'refl=\"$name_refl\"\' \'clouds=\"Observed\"\' \'site=\"$city-$subloc\"\' plots/plot_comparison_shade.ncl"; } } } } if ($plot) { # Plotting - Container experiments print "## Running container experiment plots...\n"; for (my $j = 0; $j < $exp_refl_size; $j++) { my $name_refl = $names_refl[$j]; for (my $k = 0; $k < $exp_shade_size; $k++) { my $name_shade = $names_shade[$k]; if ($cflag eq 'F') { for (my $l = 0; $l < $exp_cloud_size; $l++) { my $name_cloud = $names_cloud[$l]; system "ncl \'infile1=\"output/$city-$subloc/container_output_$city-$subloc-Coffee-$name_refl-$name_shade-$name_cloud.txt\"\' \'infile2=\"output/$city-$subloc/container_output_$city-$subloc-Bucket-$name_refl-$name_shade-$name_cloud.txt\"\' \'infile3=\"output/$city-$subloc/container_output_$city-$subloc-Barrel-$name_refl-$name_shade-$name_cloud.txt\"\' \'outfile1=\"plots/output/comparison_containers/temperature-$city-${subloc}-$name_refl-$name_shade-$name_cloud-containers\"\' \'outfile2=\"plots/output/comparison_containers/waterheight-$city-${subloc}-$name_refl-$name_shade-$name_cloud-containers\"\' \'refl=\"$name_refl\"\' \'shade=\"$name_shade\"\' \'clouds=\"$name_cloud\"\' \'site=\"$city-$subloc\"\' plots/plot_comparison_container.ncl"; } } else { system "ncl \'infile1=\"output/$city-$subloc/container_output_$city-$subloc-Coffee-$name_refl-$name_shade.txt\"\' \'infile2=\"output/$city-$subloc/container_output_$city-$subloc-Bucket-$name_refl-$name_shade.txt\"\' \'infile3=\"output/$city-$subloc/container_output_$city-$subloc-Barrel-$name_refl-$name_shade.txt\"\' \'outfile1=\"plots/output/comparison_containers/temperature-$city-${subloc}-$name_refl-$name_shade-containers\"\' \'outfile2=\"plots/output/comparison_containers/waterheight-$city-${subloc}-$name_refl-$name_shade-containers\"\' \'refl=\"$name_refl\"\' \'shade=\"$name_shade\"\' \'clouds=\"Observed\"\' \'site=\"$city-$subloc\"\' plots/plot_comparison_container.ncl"; } } } if ($cflag eq 'T') { my $ncols = 12; if ($sflag eq 'F') { $ncols = 11 } # Plotting - Cloud fraction and precipitation print "## Running cloud fraction and precipitation plots...\n"; system "ncl \'infile=\"input/container_input_$city-$subloc.txt\"\' \'outfile=\"plots/output/clouds_precip/clouds-precip-$city-${subloc}\"\' \'site=\"$city-$subloc\"\' \'ncols=$ncols\' plots/plot_clouds_precip.ncl"; # Plotting - Temperature and RH print "## Running temperature and RH plots...\n"; system "ncl \'infile=\"input/container_input_$city-$subloc.txt\"\' \'outfile=\"plots/output/temp_and_rh/temp-rh-$city-${subloc}\"\' \'site=\"$city-$subloc\"\' \'ncols=$ncols\' plots/plot_temp_rh.ncl"; } else { my $ncols = 9; if ($sflag eq 'F') { $ncols = 8 } # Plotting - Temperature and RH print "## Running temperature and RH plots...\n"; system "ncl \'infile=\"input/container_input_$city-$subloc.txt\"\' \'outfile=\"plots/output/temp_and_rh/temp-rh-$city-${subloc}\"\' \'site=\"$city-$subloc\"\' \'ncols=$ncols\' plots/plot_temp_rh.ncl"; } }
}
if ($plot)
{
# Plotting - City experiments
print "# Running city experiment plots...\n";
for (my $i = 0; $i < $exp_cont_size; $i++)
{
my $name_cont = $names_cont[$i];
for (my $j = 0; $j < $exp_refl_size; $j++)
{
my $name_refl = $names_refl[$j];
for (my $k = 0; $k < $exp_shade_size; $k++)
{
my $name_shade = $names_shade[$k];
if ($cflag eq 'F')
{
for (my $l = 0; $l < $exp_cloud_size; $l++)
{
my $name_cloud = $names_cloud[$l];
system "ncl \'infile1=\"output/Veracruz-Bachilleres/container_output_Veracruz-Bachilleres-$name_cont-$name_refl-$name_shade-$name_cloud.txt\"\' \'infile2=\"output/RioBlanco-RioBlanco/container_output_RioBlanco-RioBlanco-$name_cont-$name_refl-$name_shade-$name_cloud.txt\"\' \'infile3=\"output/Puebla-Carmen/container_output_Puebla-Carmen-$name_cont-$name_refl-$name_shade-$name_cloud.txt\"\' \'outfile1=\"plots/output/comparison_cities/temperature-$name_cont-$name_refl-$name_shade-$name_cloud\"\' \'outfile2=\"plots/output/comparison_cities/waterheight-$name_cont-$name_refl-$name_shade-$name_cloud\"\' \'shade=\"$name_shade\"\' \'refl=\"$name_refl\"\' \'clouds=\"$name_cloud\"\' \'container=\"$name_cont\"\' plots/plot_comparison_city.ncl";
}
}
else
{
system "ncl \'infile1=\"output/Veracruz-Bachilleres/container_output_Veracruz-Bachilleres-$name_cont-$name_refl-$name_shade.txt\"\' \'infile2=\"output/RioBlanco-RioBlanco/container_output_RioBlanco-RioBlanco-$name_cont-$name_refl-$name_shade.txt\"\' \'infile3=\"output/Puebla-Carmen/container_output_Puebla-Carmen-$name_cont-$name_refl-$name_shade.txt\"\' \'outfile1=\"plots/output/comparison_cities/temperature-$name_cont-$name_refl-$name_shade\"\' \'outfile2=\"plots/output/comparison_cities/waterheight-$name_cont-$name_refl-$name_shade\"\' \'shade=\"$name_shade\"\' \'refl=\"$name_refl\"\' \'clouds=\"Observed\"\' \'container=\"$name_cont\"\' plots/plot_comparison_city.ncl";
}
}
}
}
}
close MTOUTFILE;
close MWOUTFILE;
Now open the summary statistics files to run some additional plots
print "# Writing summary statistics files and plots...\n";
my %data;
my @tvals = qw/ TMEAN TMAX TMIN CMEAN CMAX CMIN /;
open MTINFILE, "<stats stats_temperature.txt"="" or="" die="" "Cannot="" open="" input="" stats="" temperature="" file:="" $!";="" open="" MWINFILE,="" "<stats="" stats_waterheight.txt"="" or="" die="" "Cannot="" open="" input="" stats="" water="" height="" file:="" $!";="" my="" $junk="<MTINFILE">;
$junk = <MWINFILE>;
while (my $tline = <MTINFILE>)
{
my $wline = <MWINFILE>;
my @tparts = split /\s+/, $tline;
my @wparts = split /\s+/, $wline;
if ($cflag eq 'F')
{
$data{$tparts[0]}->{$tparts[1]}->{$tparts[2]}->{$tparts[3]}->{$tparts[4]}->{'TMEAN'} = $tparts[5];
$data{$tparts[0]}->{$tparts[1]}->{$tparts[2]}->{$tparts[3]}->{$tparts[4]}->{'TMAX'} = $tparts[6];
$data{$tparts[0]}->{$tparts[1]}->{$tparts[2]}->{$tparts[3]}->{$tparts[4]}->{'TMIN'} = $tparts[7];
$data{$tparts[0]}->{$tparts[1]}->{$tparts[2]}->{$tparts[3]}->{$tparts[4]}->{'CMEAN'} = $tparts[8];
$data{$tparts[0]}->{$tparts[1]}->{$tparts[2]}->{$tparts[3]}->{$tparts[4]}->{'CMAX'} = $tparts[9];
$data{$tparts[0]}->{$tparts[1]}->{$tparts[2]}->{$tparts[3]}->{$tparts[4]}->{'CMIN'} = $tparts[10];
$data{$wparts[0]}->{$wparts[1]}->{$wparts[2]}->{$wparts[3]}->{$tparts[4]}->{'WH'} = $wparts[5];
}
else
{
$data{$tparts[0]}->{$tparts[1]}->{$tparts[2]}->{$tparts[3]}->{'TMEAN'} = $tparts[4];
$data{$tparts[0]}->{$tparts[1]}->{$tparts[2]}->{$tparts[3]}->{'TMAX'} = $tparts[5];
$data{$tparts[0]}->{$tparts[1]}->{$tparts[2]}->{$tparts[3]}->{'TMIN'} = $tparts[6];
$data{$tparts[0]}->{$tparts[1]}->{$tparts[2]}->{$tparts[3]}->{'CMEAN'} = $tparts[7];
$data{$tparts[0]}->{$tparts[1]}->{$tparts[2]}->{$tparts[3]}->{'CMAX'} = $tparts[8];
$data{$tparts[0]}->{$tparts[1]}->{$tparts[2]}->{$tparts[3]}->{'CMIN'} = $tparts[9];
$data{$wparts[0]}->{$wparts[1]}->{$wparts[2]}->{$wparts[3]}->{'WH'} = $wparts[4];
}
}
close MTINFILE;
close MWINFILE;
foreach my $name_refl (@names_refl)
{
foreach my $name_shade (@names_shade)
{
open WTEMPFILE, ">temp_plot_wh.txt" or die "Cannot open temp plot file: $!";
foreach my $tval (@tvals)
{
open TTEMPFILE, ">temp_plot_t_$tval.txt" or die "Cannot open temp plot file: $!";
foreach my $name_cont (@names_cont)
{
if ($cflag eq 'F')
{
foreach my $name_cloud (@names_cloud)
{
printf TTEMPFILE "%6.2f %6.2f %6.2f ", $data{"$cities[0]-$sublocs[0]"}->{$name_cont}->{$name_refl}->{$name_shade}->{$name_cloud}->{$tval}, $data{"$cities[1]-$sublocs[1]"}->{$name_cont}->{$name_refl}->{$name_shade}->{$name_cloud}->{$tval}, $data{"$cities[2]-$sublocs[2]"}->{$name_cont}->{$name_refl}->{$name_shade}->{$name_cloud}->{$tval};
}
}
else
{
printf TTEMPFILE "%6.2f %6.2f %6.2f ", $data{"$cities[0]-$sublocs[0]"}->{$name_cont}->{$name_refl}->{$name_shade}->{$tval}, $data{"$cities[1]-$sublocs[1]"}->{$name_cont}->{$name_refl}->{$name_shade}->{$tval}, $data{"$cities[2]-$sublocs[2]"}->{$name_cont}->{$name_refl}->{$name_shade}->{$tval};
}
}
print TTEMPFILE "\n";
close TTEMPFILE;
# Plotting - Temperature bar plots if ($plot) { if ($cflag eq 'F') { foreach my $name_cloud (@names_cloud) { system "ncl \'infile=\"temp_plot_t_$tval.txt\"\' \'outfile=\"plots/output/bar_temperature/$tval-$name_refl-$name_shade-$name_cloud\"\' \'shade=\"$name_shade\"\' \'refl=\"$name_refl\"\' \'clouds=\"$name_cloud\"\' plots/plot_bar_t.ncl"; } } else { system "ncl \'infile=\"temp_plot_t_$tval.txt\"\' \'outfile=\"plots/output/bar_temperature/$tval-$name_refl-$name_shade\"\' \'shade=\"$name_shade\"\' \'refl=\"$name_refl\"\' \'clouds=\"Observed\"\' plots/plot_bar_t.ncl"; } } } # Plotting - Summary DTR if ($plot) { if ($cflag eq 'F') { foreach my $name_cloud (@names_cloud) { system "ncl \'infile1=\"temp_plot_t_TMAX.txt\"\' \'infile2=\"temp_plot_t_TMIN.txt\"\' \'outfile=\"plots/output/bar_dtr/DTR-$name_refl-$name_shade-$name_cloud\"\' \'shade=\"$name_shade\"\' \'refl=\"$name_refl\"\' \'clouds=\"$name_cloud\"\' plots/plot_bar_dtr.ncl"; } } else { system "ncl \'infile1=\"temp_plot_t_TMAX.txt\"\' \'infile2=\"temp_plot_t_TMIN.txt\"\' \'outfile=\"plots/output/bar_dtr/DTR-$name_refl-$name_shade\"\' \'shade=\"$name_shade\"\' \'refl=\"$name_refl\"\' \'clouds=\"Observed\"\' plots/plot_bar_dtr.ncl"; } } foreach my $name_cont (@names_cont) { if ($cflag eq 'F') { foreach my $name_cloud (@names_cloud) { printf WTEMPFILE "%6.2f %6.2f %6.2f ", $data{"$cities[0]-$sublocs[0]"}->{$name_cont}->{$name_refl}->{$name_shade}->{$name_cloud}->{'WH'}, $data{"$cities[1]-$sublocs[1]"}->{$name_cont}->{$name_refl}->{$name_shade}->{$name_cloud}->{'WH'}, $data{"$cities[2]-$sublocs[2]"}->{$name_cont}->{$name_refl}->{$name_shade}->{$name_cloud}->{'WH'}; } } else { printf WTEMPFILE "%6.2f %6.2f %6.2f ", $data{"$cities[0]-$sublocs[0]"}->{$name_cont}->{$name_refl}->{$name_shade}->{'WH'}, $data{"$cities[1]-$sublocs[1]"}->{$name_cont}->{$name_refl}->{$name_shade}->{'WH'}, $data{"$cities[2]-$sublocs[2]"}->{$name_cont}->{$name_refl}->{$name_shade}->{'WH'}; } # Plotting - Summary frequency distribution of temperature by city if ($plot) { if ($cflag eq 'F') { foreach my $name_cloud (@names_cloud) { system "ncl \'infile1=\"output/$cities[0]-$sublocs[0]/stats_$cities[0]-$sublocs[0]-$name_cont-$name_refl-$name_shade-$name_cloud.txt\"\' \'infile2=\"output/$cities[1]-$sublocs[1]/stats_$cities[1]-$sublocs[1]-$name_cont-$name_refl-$name_shade-$name_cloud.txt\"\' \'infile3=\"output/$cities[2]-$sublocs[2]/stats_$cities[2]-$sublocs[2]-$name_cont-$name_refl-$name_shade-$name_cloud.txt\"\' \'outfile=\"plots/output/freq_distribution/tmean-$name_cont-$name_refl-$name_shade-$name_cloud\"\' \'shade=\"$name_shade\"\' \'refl=\"$name_refl\"\' \'clouds=\"$name_cloud\"\' \'container=\"$name_cont\"\' plots/plot_freq_dist_t.ncl"; } } else { system "ncl \'infile1=\"output/$cities[0]-$sublocs[0]/stats_$cities[0]-$sublocs[0]-$name_cont-$name_refl-$name_shade.txt\"\' \'infile2=\"output/$cities[1]-$sublocs[1]/stats_$cities[1]-$sublocs[1]-$name_cont-$name_refl-$name_shade.txt\"\' \'infile3=\"output/$cities[2]-$sublocs[2]/stats_$cities[2]-$sublocs[2]-$name_cont-$name_refl-$name_shade.txt\"\' \'outfile=\"plots/output/freq_distribution/tmean-$name_cont-$name_refl-$name_shade\"\' \'shade=\"$name_shade\"\' \'refl=\"$name_refl\"\' \'clouds=\"Observed\"\' \'container=\"$name_cont\"\' plots/plot_freq_dist_t.ncl"; } } } print WTEMPFILE "\n"; close WTEMPFILE; # Plotting - Summary water height distribution if ($plot) { if ($cflag eq 'F') { foreach my $name_cloud (@names_cloud) { system "ncl \'infile=\"temp_plot_wh.txt\"\' \'outfile=\"plots/output/bar_waterheight/wh-$name_refl-$name_shade-$name_cloud\"\' \'shade=\"$name_shade\"\' \'refl=\"$name_refl\"\' \'clouds=\"$name_cloud\"\' plots/plot_bar_wh.ncl"; } } else { system "ncl \'infile=\"temp_plot_wh.txt\"\' \'outfile=\"plots/output/bar_waterheight/wh-$name_refl-$name_shade\"\' \'shade=\"$name_shade\"\' \'refl=\"$name_refl\"\' \'clouds=\"Observed\"\' plots/plot_bar_wh.ncl"; } } }
}
Clean up
print "# Cleaning up...\n";
system "rm temp_plot_t_*.txt temp_plot_wh.txt param_file.txt";
print "# Finished!\n";
Subroutines
sub daily_energy_balance_stats
{
my ($einfile, $eoutfile, $soutfile, $citysubloc, $name_cont, $name_refl, $name_shade, $cflag, $name_cloud) = @_;
# Construct average daily energy balance and some other statistics
my %ebal;
my %tw;
my %tc;
my %ta;
my %evap;
my $whc = 0;
open EINFILE, "<$einfile" or die "Cannot open input file for energy balance: $!";
for (my $x = 0; $x < 7; $x++)
{
my $junk = <EINFILE>; # Header lines
}
while (my $line = <EINFILE>)
{
chomp ($line);
my @parts = split /\s+/, $line;
my $md = $parts[1].$parts[2];
my $hour = $parts[3];
if ($parts[24] != $miss)
{
push (@{$ebal{$hour}->{'SWDOWN'}}, $parts[4]);
push (@{$ebal{$hour}->{'SWSIDE'}}, $parts[5]);
push (@{$ebal{$hour}->{'LWDOWN'}}, $parts[6]);
push (@{$ebal{$hour}->{'LWUP'}}, $parts[7]);
push (@{$ebal{$hour}->{'LWIN'}}, $parts[8]);
push (@{$ebal{$hour}->{'LWOUT'}}, $parts[9]);
push (@{$ebal{$hour}->{'SHUP'}}, $parts[10]);
push (@{$ebal{$hour}->{'SHOUT'}}, $parts[11]);
push (@{$ebal{$hour}->{'LHUP'}}, $parts[12]);
push (@{$ebal{$hour}->{'GHW'}}, $parts[13]);
push (@{$ebal{$hour}->{'GHC'}}, $parts[14]);
push (@{$ebal{$hour}->{'COND'}}, $parts[15]);
push (@{$ebal{$hour}->{'BAL_W'}}, $parts[16]);
push (@{$ebal{$hour}->{'BAL_C'}}, $parts[17]);
my $sw = $parts[4] + $parts[5];
push (@{$ebal{$hour}->{'SW'}}, $sw);
my $lw = $parts[6] - $parts[7] + $parts[8] - $parts[9];
push (@{$ebal{$hour}->{'LW'}}, $lw);
my $sh = -1.*($parts[10] + $parts[11]);
push (@{$ebal{$hour}->{'SH'}}, $sh);
push (@{$evap{$md}}, $parts[19]);
push (@{$ebal{$hour}->{'TG'}}, $parts[21]);
push (@{$ebal{$hour}->{'TA'}}, $parts[22]);
push (@{$ebal{$hour}->{'TCON'}}, $parts[23]);
push (@{$ebal{$hour}->{'TW'}}, $parts[24]);
push (@{$tw{$md}}, $parts[24]);
push (@{$tc{$md}}, $parts[23]);
push (@{$ta{$md}}, $parts[22]);
if ($parts[20] >= 0.015)
{
$whc++;
}
}
}
close EINFILE;
# Write out daily average energy balance file open EOUTFILE, ">$eoutfile" or die "Cannot open output file for energy balance: $!"; print EOUTFILE "Hour,SWDOWN,SWSIDE,LWDOWN,LWUP,LWIN,LWOUT,SHUP,SHOUT,LHUP,GHW,GHC,COND,BAL_W,BAL_C,SW,LW,SH,TG,TA,TCON,TW\n"; for (my $hr = 0; $hr < 24; $hr++) { printf EOUTFILE "%02s ", $hr; my @terms = qw/ SWDOWN SWSIDE LWDOWN LWUP LWIN LWOUT SHUP SHOUT LHUP GHW GHC COND BAL_W BAL_C SW LW SH TG TA TCON TW /; foreach my $term (@terms) { my $estats = Statistics::Descriptive::Full->new(); $estats->add_data(@{$ebal{$hr}->{$term}}); my $mean = $estats->mean(); printf EOUTFILE "%8.2f ", $mean; } print EOUTFILE "\n"; } close EOUTFILE; # Calculate some basic statistics and output to file open SOUTFILE, ">$soutfile" or die "Cannot open output file for stats: $!"; my %tastats; $tastats{'TMEAN'} = Statistics::Descriptive::Full->new(); $tastats{'TMAX'} = Statistics::Descriptive::Full->new(); $tastats{'TMIN'} = Statistics::Descriptive::Full->new(); $tastats{'CMEAN'} = Statistics::Descriptive::Full->new(); $tastats{'CMAX'} = Statistics::Descriptive::Full->new(); $tastats{'CMIN'} = Statistics::Descriptive::Full->new(); $tastats{'TAMEAN'} = Statistics::Descriptive::Full->new(); $tastats{'EVAP'} = Statistics::Descriptive::Full->new(); foreach my $md (sort keys %tw) { my $tstats = Statistics::Descriptive::Full->new(); my $astats = Statistics::Descriptive::Full->new(); my $cstats = Statistics::Descriptive::Full->new(); $tstats->add_data(@{$tw{$md}}); $astats->add_data(@{$ta{$md}}); $cstats->add_data(@{$tc{$md}}); $tastats{'EVAP'}->add_data(@{$evap{$md}}); my $tmean = $tstats->mean(); my $tmax = $tstats->max(); my $tmin = $tstats->min(); my $tamean = $astats->mean(); my $cmean = $cstats->mean(); my $cmax = $cstats->max(); my $cmin = $cstats->min(); if (scalar(@{$tw{$md}}) == 24) { $tastats{'TMEAN'}->add_data($tmean); $tastats{'TMAX'}->add_data($tmax); $tastats{'TMIN'}->add_data($tmin); $tastats{'CMEAN'}->add_data($cmean); $tastats{'CMAX'}->add_data($cmax); $tastats{'CMIN'}->add_data($cmin); $tastats{'TAMEAN'}->add_data($tamean); } } print SOUTFILE "Average Water Daily Temperature (Mean, Max, Min): \n"; my $tmean = $tastats{'TMEAN'}->mean(); my $tmax = $tastats{'TMAX'}->mean(); my $tmin = $tastats{'TMIN'}->mean(); my $cmean = $tastats{'CMEAN'}->mean(); my $cmax = $tastats{'CMAX'}->mean(); my $cmin = $tastats{'CMIN'}->mean(); my $tamean = $tastats{'TAMEAN'}->mean(); my $evapmean = $tastats{'EVAP'}->mean(); printf SOUTFILE "%6.2f %6.2f %6.2f\n", $tmean, $tmax, $tmin; print SOUTFILE "Average Container Daily Temperature (Mean, Max, Min): \n"; printf SOUTFILE "%6.2f %6.2f %6.2f\n", $cmean, $cmax, $cmin; if ($cflag eq 'F') { printf MTOUTFILE "%-30s %-10s %-15s %-10s %-15s %-6.2f %-6.2f %-6.2f %-6.2f%-6.2f %-6.2f %-6.2f %-6.4f\n", $citysubloc, $name_cont, $name_refl, $name_shade, $name_cloud, $tmean, $tmax, $tmin, $cmean, $cmax, $cmin, $tamean, $evapmean; } else { printf MTOUTFILE "%-30s %-10s %-15s %-10s %-6.2f %-6.2f %-6.2f %-6.2f %-6.2f %-6.2f %-6.2f %-6.4f\n", $citysubloc, $name_cont, $name_refl, $name_shade, $tmean, $tmax, $tmin, $cmean, $cmax, $cmin, $tamean, $evapmean; } print SOUTFILE "Number of days water level > 15 mm: \n"; $whc = $whc / 24; printf SOUTFILE "%6.2f\n", $whc; if ($cflag eq 'F') { printf MWOUTFILE "%-30s %-10s %-15s %-10s %-15s %-6.2f\n", $citysubloc, $name_cont, $name_refl, $name_shade, $name_cloud, $whc; } else { printf MWOUTFILE "%-30s %-10s %-15s %-10s %-6.2f\n", $citysubloc, $name_cont, $name_refl, $name_shade, $whc; } print SOUTFILE "Frequency Distribution: \n"; my @bins = (10,12,14,16,18,20,22,24,26,28,30,32,34,36,38,40); my $fd = $tastats{'TMEAN'}->frequency_distribution_ref(\@bins); for (sort {$a <=> $b} keys %{$fd}) { print SOUTFILE "$_ $fd->{$_}\n"; } close SOUTFILE;
}
