#EFED PRZM/EXAMS Shell version 4 #original developer: Ian Kennedy # Revisions for EFED PRZM/EXAMS Shell version 4: # EFED PRZM/EXAMS Shell version 4 (file name PE4.pl dated January 8, 2003) # revised from EFED PRZM/EXAMS Shell version 1.2 # Revisions made by Dirk Young and Marietta Echeverria: # 1. Forced farm pond throughflow to zero by setting EVAP, NPSFLO, NPSED, STFLO, RAIN to zero. # 2. Revised code to read instantaneous peak concentration (rather than 24 hr average). # This was necessary because EXAMS 2.98.04 changed the location of this info. # 3. Revised code to allow EXAMS 2.98.04 to read metfile. # 4. Removed average flow option from window interface since it is no longer used. But then # put it back because of bomb throwers. # 5. Corrected heading in pzr file to read as ppm (instead of ppb). # 6. do not delete fgetsexp.xms # Revisions of 8Aug2003: # pulled the writes out of saveRun and put them in a new sub called # writeChemRunInputs(), and had this called by saveRun. Also call # this new routine from table20 so a copy of the .pzr file is appended # to the output file. use Tk; use Tk::DialogBox; use Tk::NoteBook; use Tk::LabEntry; use Tk::Checkbutton; use Tk::Listbox; #use Tk::Radiobutton; use Data::Dumper; #use strict; #$^W++; # Turn on warnings ##### Global Configuration $datemessage = "pe4.pl - 8-August-2003"; # Paths -- This is the location of the input files (path1 or path 2) and # of the PRZM and EXAMS executables. May be modified. # -- if path1 is not available it will use path2 for the input files. # The working directory is set when run files are opened or saved. my $path2 = 'F:\user\share\models\inputs'; my $path1 = "C:\\models\\INPUTS"; my $PRZMcmd = 'C:\MODELS\PRZM3\PRZM312.exe'; my $EXAMScmd = 'C:\MODELS\EXAMS\EXAMS'; my $filepath = $path1; if(opendir DIR, $filepath) { closedir DIR; } else { $filepath = $path2; } my $oldOS = 0; my @errors; # error messages (mostly where there were once "die" statements) my %C = (); my %P=(); my %Hydrol=(); my @ex; my @AppInt; my (@Tmp_mnth, @Wnd_mnth); # Average temperatures and windspeeds from metfile my ($file, $metfilename, $PRZMenv, $EXAMSenv, $path) = ("","","","","C:\\"); my ($startdate, $enddate); my ($benthic); # flags used by checkboxes my ($nApps, $firstYr, $lastYr, $outputname); # A list of variables in the main input file which need to transferred to the # %P hash before writing the PRZM *.inp file: my @CPRZMvars = ("CAM", "DEPI", "TAPP", "APPEFF", "FILTRA", "IPSCND","DRFT", "UPTKF", "PLVKRT", "PLDKRT", "FEXTRC", "PSTNAM", "DAIR", "ENPY"); # note that the %P variables KD and HENRYK are also set from %C, but these # have different names in %C (KPS and HENRY) and so are set separately. $C{FEXTRC} = 0.5; # this is a default value. #### my $DEBUG = 1; # turn off (0) if you don't want to see internals my $mainwindow; OS: { # Check if we're running WinNT. Win9x seems to not handle IPC # adequately. This sets a global variable "oldOS" which is used # by the subs that run PRZM and EXAMS. This section can be removed # when using a non-MS operating system and leaving oldOS set to zero # should work. my @OSversion = &Win32::GetOSVersion; print $OSversion[4], "\n"; $oldOS = 1 unless ($OSversion[4] == 2); } # Begin main program MAIN: { debug("+MAIN"); $mainwindow = MainWindow->new(); my $menubar = $mainwindow->Frame()->grid(-row => 0, -col => 0, -columnspan => 3, -sticky => 'nw'); my $filemenu = $menubar->Menubutton(-text => 'File'); $filemenu->command( -label => 'Open...', -command => \&loadRun ); $filemenu->command( -label => 'Save...', -command => \&saveRun ); $filemenu->separator(); $filemenu->command(-label => 'Quit', -command => sub {exit;} ); $filemenu->pack(-side => 'left'); my $editmenu = $menubar->Menubutton(-text =>'Edit'); $editmenu->command(-label => 'More PRZM paramters...', -command => \&setAdditional ); $editmenu->command(-label => 'Hydrolysis...', -command => \&setHydrol ); $editmenu->command(-label=>'PRZM env...', -command => \&editPRZMvars); $editmenu->command(-label=>'EXAMS commands...', -command => \&editEXAMScommands); $editmenu->pack(-side => 'left'); my $left = $mainwindow->Frame->grid(-row => 1, -col => 0,-sticky => 'nw'); my $right = $mainwindow->Frame->grid(-row => 1, -col => 1, -sticky => 'nw'); my $bottom = $mainwindow->Frame->grid(-row => 2, -col => 0, -columnspan => 3,-sticky => 'nw'); ### Left panel my($metlabel, $PRZlabel, $EXMlabel) = ("Metfile", "PRZM scenario", "EXAMS environment"); my @metfilelist = GetMetFiles(); my $metFrame = $left->Frame->pack(-side => "top", -anchor => "nw", -pady => 1); $metFrame->Optionmenu( -options => \@metfilelist, -width=>16, -textvariable => \$metfilename, -command => sub{($startdate, $enddate) = findmetbounds(); $firstYr = substr($startdate,-2,2); $lastYr = substr($enddate,-2,2);} )->pack(-side => "left", -anchor => "nw"); $metFrame->Label(-textvariable=>\$metlabel)->pack(-side => "left", -anchor => "nw", -fill=>'y'); $metFrame->Button(text => "...from scenario", -command => \&FindMetfile )->pack(-side => "right", -anchor => "nw"); my @PRZMfilelist = GetPRZMenvs(); $PRZMenv = $PRZMfilelist[5]; # default value my $PRZMframe = $left->Frame->pack(-side => "top", -anchor => "nw", -pady => 1); $PRZMframe->Optionmenu( -options => \@PRZMfilelist, -width=>16, -textvariable => \$PRZMenv, -command => sub{my $x = $PRZMenv;ReadPRZMenv($x);}, )->pack(-side => "left", -anchor => "nw"); $PRZMframe->Label(-textvariable=>\$PRZlabel)->pack(-side => "left", -anchor => "nw", -fill=>'y'); my @EXAMSfilelist = GetEXAMSenvs(); my $EXMframe=$left->Frame->pack(-side => "top", -anchor => "nw", -pady => 1); $EXMframe->Optionmenu( -options => \@EXAMSfilelist, -width=>16, -textvariable => \$EXAMSenv, )->pack(-side => "left", -anchor => "nw"); $EXMframe->Label(-textvariable=>\$EXMlabel)->pack(-side => "left", -anchor => "nw", -fill=>'y'); my $fieldFrame = $left->Frame->pack(-side => "top", -pady => 1, -anchor => "nw"); $fieldFrame->Label(-text => "Field size:" )->pack(-side => "left", -anchor => "nw", -pady => 4); # Make sure $C{IR} is set to a valid value (Default, IR or Pond) $C{IR} = "Default" unless($C{IR} =~ /IR|Pond/); $fieldFrame->Radiobutton(-text => "IR", -value => "IR", -variable => \$C{IR} )->pack(-side => "left", -anchor => "nw", -pady => 1); $fieldFrame->Radiobutton(-text=> "Pond", -value => "Pond", -variable => \$C{IR} )->pack(-side => "left", -anchor => "e", -pady => 1); $fieldFrame->Radiobutton(-text=> "Default",-value => "Default", -variable => \$C{IR} )->pack(-side => "right", -anchor => "center", -pady => 1); my $runfFrame = $left->Frame->pack(-side => "top", -pady => 1, -anchor => "nw"); $runfFrame->Label(-text => "Runoff flow:" )->pack(-side => "left", -anchor => "nw", -pady => 4); $runfFrame->Radiobutton(-text => "Monthly", -value => "monthly", -variable => \$C{RUNOFF} )->pack(-side => "left", -anchor => "nw", -pady => 1); # monthly option (control only) removed 07/01/03 IK. put back in by dfy $runfFrame->Radiobutton(-text=> "Overall (Sure?)", -value => "total", -variable => \$C{RUNOFF} )->pack(-side => "left", -anchor => "e", -pady => 1); $runfFrame->Radiobutton(-text=> "None",-value => "none", -variable => \$C{RUNOFF} )->pack(-side => "right", -anchor => "center", -pady => 1); # $left->Checkbutton(-text => "Use runoff for res. input", # -variable => \$C{RUNOFF}, # -onvalue => "monthly", -offvalue => "none" # )->pack(-side => "top", # -anchor => "nw", -pady => 1); $left->LabEntry(-label => "Output filename", -labelPack => [-side => "right", -anchor => "w"], -width => 24, -bg => 'white', -textvariable => \$outputname)->pack(-side => "top", -pady => 1, -anchor => "nw"); $left->LabEntry(-label => "Chemical name", -labelPack => [-side => "right", -anchor => "w"], -width => 16, -bg => 'white', -textvariable => \$C{PSTNAM})->pack(-side => "top", -pady => 1, -anchor => "nw"); $left->LabEntry(-label => "Molecular Weight", -labelPack => [-side => "right", -anchor => "w"], -width => 16, -bg => 'white', -validate => "all", -validatecommand => \&val8, -textvariable => \$C{MWT})->pack(-side => "top", -pady => 1, -anchor => "nw"); $left->LabEntry(-label => "Henry's Law Const. (atm m^3/mol) ", -labelPack => [-side => "right", -anchor => "w"], -width => 16, -bg => 'white', -validate => "all", -validatecommand => \&val8, -textvariable => \$C{HENRY})->pack(-side => "top", -pady => 1, -anchor => "nw"); $left->LabEntry(-label => "Vapour Pressure (torr)", -labelPack => [-side => "right", -anchor => "w"], -width => 16, -bg => 'white', -validate => "all", -validatecommand => \&val8, -textvariable => \$C{VAPR})->pack(-side => "top", -pady => 1, -anchor => "nw"); $left->LabEntry(-label => "Solubility (mg/L)", -labelPack => [-side => "right", -anchor => "w"], -width => 16, -bg => 'white', -validate => "all", -validatecommand => \&val8, -textvariable => \$C{SOL})->pack(-side => "top", -pady => 1, -anchor => "nw"); $left->LabEntry(-label => "Kd", -labelPack => [-side => "right", -anchor => "w"], -width => 16, -bg => 'white', -validate => "all", -validatecommand => \&val8, -textvariable => \$C{KD})->pack(-side => "top", -pady => 1, -anchor => "nw"); $left->LabEntry(-label => "Koc", -labelPack => [-side => "right", -anchor => "w"], -width => 16, -bg => 'white', -validate => "all", -validatecommand => \&val8, -textvariable => \$C{KOC})->pack(-side => "top", -pady => 1, -anchor => "nw"); # Right Panel my $CAMframe = $right->Frame->pack(-side => "top", -anchor => "nw", -pady => 1); my @CAMlist = ("1","2","4","5","6","7","8"); my $CAMlabel = "CAM"; $CAMframe->Optionmenu( -options => \@CAMlist, -width=>6, -textvariable => \$C{CAM}, )->pack(-side => "left", -anchor => "nw"); $CAMframe->Label(-textvariable=>\$CAMlabel)->pack(-side => "left", -anchor => "nw", -fill=>'y'); $right->LabEntry(-label => "Incorp. Depth (cm)", -labelPack => [-side => "right", -anchor => "w"], -width => 16, -bg => 'white', -validate => "all", -validatecommand => \&val5, -textvariable => \$C{DEPI})->pack(-side => "top", -pady => 1, -anchor => "nw"); $right->LabEntry(-label => "App. Rate (kg a.i./ha)", -labelPack => [-side => "right", -anchor => "w"], -width => 16, -bg => 'white', -textvariable => \$C{TAPP}, -validate => "all", -validatecommand => \&val6, )->pack(-side => "top", -pady => 1, -anchor => "nw"); $right->LabEntry(-label => "App. Efficiency (fraction)", -labelPack => [-side => "right", -anchor => "w"], -width => 16, -bg => 'white', -validate => "all", -validatecommand => \&val5, -textvariable => \$C{APPEFF})->pack(-side => "top", -pady => 1, -anchor => "nw"); $right->LabEntry(-label => "Spray Drift (fraction)", -labelPack => [-side => "right", -anchor => "w"], -width => 16, -bg => 'white', -validate => "all", -validatecommand => \&val5, -textvariable => \$C{DRFT})->pack(-side => "top", -pady => 1, -anchor => "nw"); $right->LabEntry(-label => "Application Date (day-mon)", -labelPack => [-side => "right", -anchor => "w"], -width => 16, -bg => 'white', -textvariable => \$C{DATE})->pack(-side => "top", -pady => 1, -anchor => "nw"); my $AppIntframe = $right->Frame->pack(-side => "top", -pady => 1, -anchor => "nw"); $nApps = @AppInt; $AppIntframe->LabEntry(-label => "Number of Applications", -labelPack => [-side => "right", -anchor => "w"], -width => 6, -bg => 'white', -validate => "all", -validatecommand => \&valApp, -textvariable => \$nApps)->pack(-side => "left", -pady => 1, -anchor => "nw"); $AppIntframe->Button(text => "Intervals...", -command => sub{ $P{NAPS}=$nApps; setAppIntervals($P{NAPS});} )->pack(-side => "right", -anchor => "ne"); my $IPSframe = $right->Frame->pack(-side => "top", -pady => 1, -anchor => "nw"); $IPSframe->Optionmenu( -options => [1,2,3], -width=>6, -variable => \$C{IPSCND}, )->pack(-side => "left", -anchor => "nw"); my $IPSlabel = "IPSCND (Record 17)"; $IPSframe->Label(-textvariable=> \$IPSlabel)->pack(-side => "left", -anchor => "nw", -fill => 'y'); my $buttonFrame = $right->Frame->pack(-side => "top", -pady => 1, -anchor => 'nw'); $buttonFrame->Button(text => "Set Hydrolysis...", -command => \&setHydrol )->pack(-side => "left", -anchor => "nw"); $buttonFrame->Button(text => "More PRZM Parameters...", -command => \&setAdditional )->pack(-side => "right", -anchor => "ne"); $right->Checkbutton(-text => "Write Benthic pore water concentrations", -variable => \$benthic)->pack(-side => "top", -anchor => "nw", -pady => 1); $right->LabEntry(-label => "Aq. Photolysis half-life (days)", -labelPack => [-side => "right", -anchor => "w"], -width => 16, -bg => 'white', -validate => "all", -validatecommand => \&val8, -textvariable => \$C{KDP})->pack(-side => "top", -pady => 1, -anchor => "nw"); $right->LabEntry(-label => "Water half-life (days)", -labelPack => [-side => "right", -anchor => "w"], -width => 16, -bg => 'white', -validate => "all", -validatecommand => \&val8, -textvariable => \$C{KBACW})->pack(-side => "top", -pady => 1, -anchor => "nw"); $right->LabEntry(-label => "Benthic half-life(days)", -labelPack => [-side => "right", -anchor => "w"], -width => 16, -bg => 'white', -textvariable => \$C{KBACS})->pack(-side => "top", -pady => 1, -anchor => "nw"); $right->LabEntry(-label => "Soil half-life (days)", -bg=>'white', -labelPack => [-side => "right", -anchor => "w"], -width => 16, -bg => 'white', -textvariable => \$C{ASM})->pack(-side => "top", -pady => 1, -anchor => "nw"); ### Bottom panel, both input & output $bottom->Button(-text => 'Run PRZM/EXAMS', -command => \&RunModels )-> grid(qw/-row 2 -column 0 -sticky nesw/); # Protect from stray signals # You could just call MainLoop() w/o the fancy eval print "$datemessage\n"; while (1) { eval MainLoop(); # Start the event processing } # Will never get here debug ("-MAIN"); } ####################################################### # Validation routines -- Make sure entrys are not too long # and reject any obviously invalid characters sub val5 { return 0 unless $_[1] =~ /[0-9\.+-eE]/ or !$_[1]; my $maxLgth = $_[0] =~ /^\./ ? 4 : 5; return length($_[0]) > $maxLgth ? 0 : 1; } sub val6 { return 0 unless $_[1] =~ /[0-9\.+-eE]/ or !$_[1]; my $maxLgth = $_[0] =~ /^\./ ? 5 : 6; return length($_[0]) > $maxLgth ? 0 : 1; } sub val8 { return 0 unless $_[1] =~ /[0-9\.+-eE]/ or !$_[1]; my $maxLgth = $_[0] =~ /^\./ ? 7 : 8; return length($_[0]) > $maxLgth ? 0 : 1; } sub valApp { return 1 if $_[4] == 0; return 0 unless $_[1] =~ /[0-9]/ or !$_[1]; return ($_[0]<1 or $_[0]>50 or $_[0] eq "") ? 0 : 1; } ######################################################### # given a filename (with path), returns a date string sub fileModTime { my $fn = shift; my @a = stat $fn; my @lt = localtime( $a[9]); my $dow = ('Sunday', 'Monday', 'Tueday', 'Wedday', 'Thuday', 'Friday', 'Satday')[@lt[6]]; my $mnth = ('January', 'February', 'March', 'April', 'May', 'June', 'July', 'August', 'September', 'October', 'November', 'December')[@lt[4]]; my $yr = 1900 + $lt[5]; my $time = sprintf("%02d:%02d:%02d", $lt[2], $lt[1], $lt[0]); return "$dow, $lt[3] $mnth $yr at $time"; } ######################################################### # save Chemical and run information. # sub saveRun { debug("+saveRun"); # Types are listed in the dialog widget my @types = (["PRZM/EXAMS Run files", '.pzr', 'TEXT'], ["All Files", "*"] ); # Uses standard file dialog for OS my $file = $mainwindow->getSaveFile(-filetypes => \@types, -initialfile => $file, -defaultextension => '.pzr'); my @path = split(/\//, $file); # get the path out of $file and use it to pop @path; # set the current working directory. my $newcwd = join("\\", @path); use Cwd 'chdir'; chdir $newcwd; print "working dir set to ", $ENV{'PWD'},"\n"; $path = $newcwd; my $filerr; # Write out in tab-delimited format open my $OUT, ">$file" or $filerr = "couldn't open $file"; if($filerr) { $mainwindow->messageBox(-icon => 'questhead', -message => $filerr, -title => 'file error', -type => 'OK'); } writeChemRunInputs($OUT); close $OUT; debug("-saveRun"); return; } sub writeChemRunInputs { my $OUT = shift; print $OUT "Output File: $outputname\n"; print $OUT "Metfile:\t$metfilename\n"; print $OUT "PRZM scenario:\t$PRZMenv\n"; print $OUT "EXAMS environment file:\t$EXAMSenv\n"; print $OUT "Chemical Name:\t$C{PSTNAM}\n"; print $OUT "Description\tVariable Name\tValue\tUnits\tComments\n"; print $OUT "Molecular weight\tmwt\t$C{MWT}\tg/mol\n"; print $OUT "Henry's Law Const.\thenry\t$C{HENRY}\tatm-m^3/mol\n"; print $OUT "Vapor Pressure\tvapr\t$C{VAPR}\ttorr\n"; print $OUT "Solubility\tsol\t$C{SOL}\tmg/L\n"; print $OUT "Kd\tKd\t$C{KD}\tmg/L\n"; print $OUT "Koc\tKoc\t$C{KOC}\tmg/L\n"; print $OUT "Photolysis half-life\tkdp\t$C{KDP}\tdays\tHalf-life\n"; print $OUT "Aerobic Aquatic Metabolism\tkbacw\t$C{KBACW}\tdays\tHalf\life\n"; print $OUT "Anaerobic Aquatic Metabolism\tkbacs\t$C{KBACS}\tdays\tHalf\life\n"; print $OUT "Aerobic Soil Metabolism\tasm\t$C{ASM}\tdays\tHalf\life\n"; for my $key(sort {$a <=> $b} keys %Hydrol) { print $OUT "Hydrolysis:\tpH $key\t$Hydrol{$key}\tdays\tHalf-life\n"; } print $OUT "Method:\tCAM\t$C{CAM}\tinteger\tSee PRZM manual\n"; print $OUT "Incorporation Depth:\tDEPI\t$C{DEPI}\tcm\n"; print $OUT "Application Rate:\tTAPP\t$C{TAPP}\tkg/ha\n"; print $OUT "Application Efficiency:\tAPPEFF\t$C{APPEFF}\tfraction\n"; print $OUT "Spray Drift\tDRFT\t$C{DRFT}\tfraction of application rate applied to pond\n"; print $OUT "Application Date\tDate\t$C{DATE}\tdd/mm or dd/mmm or dd-mm or dd-mmm\n"; for(my $i; $i<@AppInt ; $i++) { print $OUT "Interval ",$i+1,"\tinterval\t$AppInt[$i]\tdays"; print $OUT "\tSet to 0 or delete line for single app.\n"; } print $OUT "Record 17:\tFILTRA\t$C{FILTRA}\n"; print $OUT "\tIPSCND\t$C{IPSCND}\n"; print $OUT "\tUPTKF\t$C{UPTKF}\n"; print $OUT "Record 18:\tPLVKRT\t$C{PLVKRT}\n"; print $OUT "\tPLDKRT\t$C{PLDKRT}\n"; print $OUT "\tFEXTRC\t$C{FEXTRC}\n"; print $OUT "Flag for Index Res. Run\tIR\t$C{IR}\n"; print $OUT "Flag for runoff calc.\tRUNOFF\t$C{RUNOFF}\tnone, monthly or total(average of entire run)\n"; } ############################################################ # Create a dialogue box for additional paramters. These include # The parameters in Records 17 and 18, and Record 26 sub setAdditional { debug("+setAdditional"); # if (not defined $addDial) { my $addDial = $mainwindow->DialogBox( -title => "Additional PRZM parameters", -buttons => ["OK"]); my $left = $addDial->Frame->grid(-row => 1, -col => 0, -sticky => 'nw'); my $right= $addDial->Frame->grid(-row => 1, -col => 1, -sticky => 'nw'); $left->LabEntry(-label => "DAIR - air diffusion (cm^2/day)", -labelPack => [-side => "right", -anchor => "w"], -width => 16, -bg => 'white', -textvariable => \$C{DAIR})->pack(-side => "top", -pady => 1, -anchor => "nw"); $left->LabEntry(-label => "ENPY - Vap. Enthalpy (kcal/mol)", -labelPack => [-side => "right", -anchor => "w"], -width => 16, -bg => 'white', -textvariable => \$C{ENPY})->pack(-side => "top", -pady => 1, -anchor => "nw"); $left->LabEntry(-label => "FILTRA - Filtration Parameter", -labelPack => [-side => "right", -anchor => "w"], -width => 16, -bg => 'white', -textvariable => \$C{FILTRA})->pack(-side => "top", -pady => 1, -anchor => "nw"); $left->LabEntry(-label => "UPTKF-Plant uptake factor", -labelPack => [-side => "right", -anchor => "w"], -width => 16, -bg => 'white', -textvariable => \$C{UPTKF})->pack(-side => "top", -pady => 1, -anchor => "nw"); $right->LabEntry(-label => "PLVKRT", -labelPack => [-side => "right", -anchor => "w"], -width => 16, -bg => 'white', -textvariable => \$C{PLVKRT})->pack(-side => "top", -pady => 1, -anchor => "nw"); $right->LabEntry(-label => "PLDKRT", -labelPack => [-side => "right", -anchor => "w"], -width => 16, -bg => 'white', -textvariable => \$C{PLDKRT})->pack(-side => "top", -pady => 1, -anchor => "nw"); $right->LabEntry(-label => "FEXTRC", -labelPack => [-side => "right", -anchor => "w"], -width => 16, -bg => 'white', -textvariable => \$C{FEXTRC})->pack(-side => "top", -pady => 1, -anchor => "nw"); my $result = $addDial->Show; # Execute the dialog box debug("-setAdditional"); } ############################################################ # Create a Hydrolysis dialogue sub setHydrol { debug("+setHydrol"); # Extract the %Hydrol hash into two local arrays, for pH (the keys) and rates my @hRate; my @pH = sort {$a <=> $b} keys %Hydrol; for(my $i=0 ; $i<3 ; $i++) { $hRate[$i] = $Hydrol{$pH[$i]}; } if (not defined $hDial) { $hDial = $mainwindow->DialogBox(-title => "Hydrolysis", -buttons => ["OK", "Cancel"]); my $left = $hDial->Frame->grid(-row => 1, -col => 0, -sticky => 'nw'); my $right= $hDial->Frame->grid(-row => 1, -col => 1, -sticky => 'nw'); $left->Label(-text => "pH", -width=>6)->pack(-side => "top", -anchor => "nw"); $right->Label(-text => "half-life (days)", -width=>12)-> pack(-side => "top", -anchor => "nw"); for(my $i=0 ; $i<3 ; $i++) { $pHwig[$i] = $left->Entry( -width => 6, -textvariable => \$pH[$i])-> pack(-side => "top", -anchor => "nw"); $hRateWig[$i] = $right->Entry( -width => 12, -textvariable => \$hRate[$i])-> pack(-side => "top", -anchor => "nw"); } } else { for(my $i=0 ; $i<3 ; $i++) { $pHwig[$i]->configure(-textvariable => \$pH[$i]); $hRateWig[$i]->configure(-textvariable => \$hRate[$i]); } } my $result = $hDial->Show; # Execute the dialog box if ($result =~ /OK/) { # Copy all the values back to %Hydrol. The keys may have changed, so # delete all values and start over %Hydrol = (); for my $i(0..2) { next if($pH[$i] eq ""); $Hydrol{$pH[$i]} = $hRate[$i]; } debug(Dumper $config); } debug("-setHydrol"); } sub FindMetfile { my $message = "Scenario does not suggest\nan MLRA or metfile."; my $met; my $wholepath = $filepath."\\PRZMenv\\".$PRZMenv; if(open PRZ, $wholepath) { while() { if(/Metfile:?\s?([Mm]et\d{1,3}[abc]?\.met|[Ww]\d{5}\.dvf)/i) { $met = "$1"; $message = "Try using $met"; last; } } } else { $message = "Couldn't open $wholepath"; } # Now display the message: my $r = $mainwindow->messageBox( -message => $message, -title => 'Metfile', -type => 'OKCancel'); if($r eq 'ok' and $met) { # Set the metfile my @mflist = GetMetFiles(); foreach my $f(@mflist) { if($f =~ /$met/i) { $metfilename = $f; ($startdate, $enddate) = findmetbounds(); $firstYr = substr($startdate,-2,2); $lastYr = substr($enddate,-2,2); last; } } } } ################################################### # Edit PRZM (%P) parameters sub editPRZMvars { my @pkeys = sort keys %P; # the following list is the variables in %P which are edited in the main window. my @noeditlist = ('PSTNAM', 'DRFT', 'FILTRA', 'CAM', 'TAPP', 'APPEFF', 'IPSCND', 'UPTKF', 'PLVKRT', 'PLDKRT', 'FEXTRC'); for(my $k=0 ; $k<@pkeys ; $k++) { foreach my $mvar(@noeditlist) { if($pkeys[$k] eq $mvar) { splice(@pkeys, $k--, 1); last; } } } my $pDial = $mainwindow->DialogBox(-title => "Edit PRZM parameters", -buttons => ["Done"]); my $left = $pDial->Frame->grid(-row => 1, -col => 0, -sticky => 'nw'); my $right = $pDial->Frame->grid(-row => 1, -col => 1, -sticky => 'nw'); # $left->Label(-text => "Name", -width=>6)->pack(-side => "top", -anchor => "nw"); # $right->Label(-text => "Value", -width=>12); my $pNameBox = $left->Listbox( -width => 10, -height => 8)-> pack(-side => "top", -anchor => "nw"); foreach (@pkeys) { $pNameBox->insert('end', $_); } my $scroll = $left->Scrollbar( -command => ['yview', $pNameBox]); $pNameBox->configure(-yscrollcommand => ['set', $scroll]); $pNameBox->pack(-side => 'left', -fill => 'both', -expand => 'yes'); $scroll->pack(-side => 'right', -fill => 'y'); my $pValWig = $right->Entry( -width => 16, -textvariable => "")-> pack(-side => "top", -anchor => "nw"); my $s; # current selection in the list box my $changebutton = $right->Button(-text=>, "Store", -command => sub{$P{$pkeys[$s]} = $pValWig->get; print "set to $P{$pkeys[$s]}\n";} )-> pack(-side => "right", -anchor => "nw"); my $editbutton = $right->Button(-text=>, "Edit", -command => sub{($s) = $pNameBox->curselection; $pValWig->configure(-textvariable=>\$P{$pkeys[$s]}); })-> pack(-side => "left", -anchor => "nw"); my $result = $pDial->Show; # Execute the dialog box } ####################### # Edit EXAMS commands # ####################### sub editEXAMScommands { my $exDial = $mainwindow->DialogBox(-title => "Edit EXAMS commands", -buttons => ["Done"]); my $top = $exDial->Frame->grid(-row => 1, -col => 0, -sticky => 'nw'); my $bottom = $exDial->Frame->grid(-row => 2, -col => 0, -sticky => 'nw'); my $exCmdBox = $top->Listbox( -width => 72, -height => 8)-> pack(-side => "top", -anchor => "nw"); foreach (@ex) { $exCmdBox->insert(0, $_); } my $scroll = $top->Scrollbar( -command => ['yview', $exCmdBox]); $exCmdBox->configure(-yscrollcommand => ['set', $scroll]); $exCmdBox->pack(-side => 'left', -fill => 'both', -expand => 'yes'); $scroll->pack(-side => 'right', -fill => 'y'); my $exEntry = $bottom->Entry( -width => 60)-> pack(-side => "left", -anchor => "nw"); my $addbutton = $bottom->Button(-text=>, "Add", -command => sub{my $x = $exEntry->get; splice(@ex,$exCmdBox->curselection,0,$x); $exCmdBox->insert('end', $x); $exEntry->delete(0,'end');} )->pack(-side => "right", -anchor => "nw"); my $deletebutton = $bottom->Button(-text=>, "Delete", -command => sub{my $x = $exCmdBox->curselection; splice(@ex,$x,1); $exCmdBox->delete($x);} )->pack(-side => "left", -anchor => "nw"); my $result = $exDial->Show; # Execute the dialog box } ################################################################ # Let the user choose more than four applications sub setAppNumber { my $appDial = $mainwindow->DialogBox(-title => "Set application number", -buttons => ["Done"]); my $number = $P{NAPS}; $appDial->LabEntry(-label => "Number of Applications", -bg=>'white', -labelPack => [-side => "right", -anchor => "w"], -width => 16, -bg => 'white', -textvariable => \$number,)->pack(-side => "top", -pady => 1, -anchor => "nw"); my $result = $appDial->Show; $result = int $result; $result = 1 if $result < 1; $result = 30 if result > 30; $P{NAPS} = $result; setAppIntervals($result) if($result > 1); debug("-setAppNumber"); return $result; } ################################################################# # Dialog box for editing applicaton intervals sub setAppIntervals { my $n = shift; # The number of applications my $appDial = $mainwindow->DialogBox(-title => "Edit application intervals", -buttons => ["Done"]); my $title = "Enter application\nintervals (days)"; $appDial->Label(-textvariable=> \$title)->pack(-side => "top", -anchor => "nw", -fill => 'y'); if($n > 1) { for my $i(1..$n-1) { $appDial->LabEntry(-label => "Interval $i:", -labelPack => [-side => "left", -anchor => "w"], -width => 6, -bg => 'white', -textvariable => \$AppInt[$i-1])->pack(-side => "top",-padx => 1, -pady=>1, -anchor => "nw"); } my $result = $appDial->Show; # Execute the dialog box } for(my $i = $#AppInt ; $i>$n-2 ; $i--) { pop @AppInt; # remove extra values from @AppInt } debug("-setAppIntervals"); } ############################################################## # Read in a run filename and set the working directory to the # directory containing it. # sub loadRun { debug("+loadRun"); # Types are listed in the dialog widget my @types = (["PRZM/EXAMS run files", '.pzr', 'TEXT'], ["Text files", '.txt', 'TEXT'], ["All Files", "*"] ); my $file = $mainwindow->getOpenFile(-filetypes => \@types); ReadFiles($file); my @path = split(/\//, $file); $file = pop @path; my $newcwd = join("\\", @path); use Cwd 'chdir'; chdir $newcwd; $path = $newcwd; debug("-loadRun"); } ################################################################### # Display the @errors array in a dialog box and then clear the array sub displayWarnings { my $message = "@errors"; my $r = $mainwindow->messageBox( -message => $message, -title => 'Warnings', -type => 'OKCancel'); @errors = (); # erase the displayed messages return $r; } ###################################################################### # three subs to get lists of filenames. One for each environment file # and one for the metfile. The metfile sub list both .awd and .met # files and removes any .awd file with the same name as a .met file # from the list. # sub GetMetFiles { opendir DIR, $filepath."\\metfiles"; my @files = grep {/\.met$|\.awd$|\.dvf$/i} readdir(DIR); foreach(@files) { tr/[A-Z]/[a-z]/; # make every name lower case this would be a problem on } # a case-sensitive system @files = sort @files; # After sorting, go through the files and eliminate any .awd files for which # there is a .met file with the same name. my ($curname, $curext) = ("",""); my ($lastname, $lastext) = split(/\./, @files[0]); for(my $i=1 ; $i<@files ; $i++) { # remove .AWD files if there is a .met ($curname, $curext) = split(/\./, @files[$i]); if($curname eq $lastname) { if($curext =~ /met/i and $lastext =~ /awd/i) { splice(@files, --$i, 1); # decrement i because indexes will shift down one } } ($lastname, $lastext) = ($curname, $curext); } return @files; } sub GetPRZMenvs { opendir DIR, $filepath."\\PRZMenv"; my @files = grep {/\.txt$/i} readdir(DIR); return sort @files; } sub GetEXAMSenvs { opendir DIR, $filepath."\\EXAMSenv"; my @files = grep {/\.exv$/i} readdir(DIR); return @files; } ###################################################################### # Print debugging information to stdout sub debug { my @msg = shift; print @msg, "\n" if $DEBUG; } ###################################################################### # Run The Models -- Called when the "Run PRZM/EXAMS button is pressed sub RunModels { print "IR = $C{IR}; Runoff = $C{RUNOFF}\n"; my $t1 = time(); SetChemParameters(); $apps = scalar(split(/\t/, $P{APD})); $P{NAPS} = $apps * ($lastYr - $firstYr + 1);# overwrite this value because # we're using the metfile data and header # info to calculate application dates if(TestData()) { # do a few checks on the data to flag some possible problems warn "$startdate - $enddate\n"; WriteRunFile(); WritePRZMinput(); if(@errors) { my $stoprun= displayWarnings(); print "stoprun = $stoprun\n"; return if($stoprun eq 'cancel'); } RunPRZM(); # see below WriteEXAfile(); # The EXAMS instructions file RunEXAMS(); # see below RunPRZM Table20("","Direct", "Water"); Table20("ben", "Direct", "Benthic") if($benthic); TimeSeries(); # Reformat and rename the XGETSEXP.XMS file print "Ran PRZM and EXAMS. Results are in $outputname.out\n"; system "del P2E-C*"; system 'del restart.prz'; system 'del p2e.exa'; # the (invalid) PRZM generated .exa file print "Total runtime ",time()-$t1, " seconds\n"; } } sub RunPRZM { my ($y,$firsty, $X) = (0,0,0); if($oldOS) { $X = system($PRZMcmd); die("coudn't run PRZM - $?\n") if($X!=0); # Error check } else { open(PRZM, $PRZMcmd.' |') or die "can't run przm\n"; print "PRZM: "; while() { if(/Simulating from /) { my $cy = substr($',7,2); $firsty = $cy if($firsty == 0); if($y != $cy) { print "$cy "; $y = $cy; } #line break every 20 years print "\n " if(($y-$firsty) % 20 == 0 and ($y-$firsty > 0) and m/Jan/); } } print "\n"; close PRZM; print "Problem running PRZM\n" if $?; } } sub RunEXAMS { my ($y, $firsty,$X) = (0,0,0); if($oldOS) { $X = system($EXAMScmd.' < pz2ex.exa'); die("coudn't run EXAMS - $?\n") if($X!=0); # Error check } else { open(EXAMS, 'c:\models\exams\exams < pz2ex.exa |') or die "can't run przm\n"; print "EXAMS: "; while() { if(/Processing com.*(\d{4})/) { my $cy = $1; $firsty = $cy if($firsty == 0); if($y != $cy) { print "$cy "; $y = $cy; } #line break every 10 years print "\n " if(($y-$firsty+1) % 10 == 0 and ($y-$firsty > 0)); } } print "\n"; close EXAMS; WriteAnex(); #system('C:\models\Anex\command_line_anex.exe'); #print "finished running Anex\n"; } } ##################################################### sub ReadFiles { my $fname = shift @_; @AppInt = (); # remove values from %C. We don't want to delete the entries, because Tk # needs them as variables. foreach my $k(keys %C) { # erase any existing data just to be safe. $C{$k} = "" ; } # These may be overwritten, but defaults are better than a previous setting ($C{IR}, $C{FEXTRC}) = ("Default", "0.5"); open(CHEM, $fname) or push @errors, "Couldn't open $fname\n"; # Read the header. The last line of the header is assumed to contain the first # instance of the word "Record" in the tab-delimited spreadsheet file. while() { chomp; chomp; s/\s+$//; # Remove any trailing whitespace if(/metfile:\s*/i) {$metfilename = $';} if(/EXAMS.*:\s*/i) {$EXAMSenv = $';} if(/PRZM.*:\s*/i) {$PRZMenv = $';} $outputname = $' if (/Output.*:\s*/i); if(/Chemical.*:\s*/i) { ($C{PSTNAM}) = split(/\t/,$'); $outputname = $C{PSTNAM} if length($outputname)<1; } last if(/Description/i); } CheckFiles(); # The three filenames are global, so no need to pass parameters ($startdate, $enddate) = findmetbounds(); $firstYr = substr($startdate,-2,2); # found when writing out "przm3.run" $lastYr = substr($enddate,-2,2); ReadChem(); # Continue reading beyond the header of the file close CHEM; ReadPRZMenv($PRZMenv); print "PRZM environment file is $PRZMenv\n"; debug("-ReadFiles -- CAM = $C{CAM}\n"); } ###################################################### # This routine checks to be sure the PRZM environment file, the EXAMS # environment file and the metfile really exist in the proper directories sub CheckFiles { my $fn; my @badfiles = ("PRZM scenario"); my @filelist = GetPRZMenvs(); foreach $fn(@filelist) { if($fn =~ /$PRZMenv/i) { pop @badfiles; last; } } # Reset PRZMenv if the file wasn't found. $PRZMenv = $filelist[0] if($badfiles[$#badfiles] =~ /PRZM/); push @badfiles, "EXAMS environment"; my @filelist = GetEXAMSenvs(); foreach $fn(@filelist) { if($fn =~ /$EXAMSenv/i) { pop @badfiles; last; } } $EXAMSenv = $filelist[0] if($badfiles[$#badfiles] =~ /EXAMS/); push @badfiles, "Metfile"; my @filelist = GetMetFiles(); foreach $fn(@filelist) { if($fn =~ /$metfilename/i) { pop @badfiles; last; } } $metfilename = $filelist[0] if($badfiles[$#badfiles] =~ /Met/); if(@badfiles > 0) { # make an alert box if there is at least one bad filename require Tk::Dialog; my $plural = @badfiles > 1 ? "s" : ""; my $temp = $"; $" = "\n"; my $message = "The file$plural:\n@badfiles\ncannot be found.\n"; $message .= "Substituting default file$plural."; my $alert = $mainwindow->Dialog(-title=>"Missing file error", -text => $message, -buttons => ["OK"]); my $result = $alert->Show; $" = $temp; } return; } ##################################################### sub ReadPRZMenv { my $name = shift @_; return unless($name); foreach my $k(keys %P) { delete $P{$k}; # $P{$k} = ""; # set all data to the empty string. this ensures all } # variables are blank to Tk as well as the rest of the script $name = $filepath."\\PRZMenv\\".$name; # add in the path open(SHEET, $name) or push @errors, "Couldn't open $name\n"; # Loop through the spreadsheet file and read all the data into a hash. # This looks for the variable name in column 2 ($fields[1]) and the data value # in column 3 ($fields[2]). If there is no variable name, or the name has # been used before (tested by whether or not there's a value for it in the # hash table) then the data value is appended to the last used name. This # does not check to see if there is a name for the first data, a blank in the # first record will cause an error. $P{WINDAY} = 0; # WINDAY is in record 16, but is not currently in the Chemical input # file or the Tk front end. Creating it here allows it to exist # and be edited as a PRZM parameter. while() { last if(/Record/i); } while() { # first make sure we have a variable name in $fields[1], and data in $fields[2]: chomp; chomp; next if(!$_); # Skip blank lines my @fields = split /\t/; if($fields[0] =~ /EXAMS/i) { $fields[1] =~ s/"//g; # remove any " symbols, EXCEL seems to insert them push(@ex, $fields[1]); # when there is a "*" in the field. } next if($fields[1] =~ /[a-z]{2}/); #Skip any line with 2 lowercase letters in the var name next if(length($fields[1])==0 and length($fields[2])==0); # skip any line w/o data if($fields[1] =~ /\s+/){ # drop anything after whitespace in the var name $fields[1] = $`; } if($fields[1] !~ /PSTNAM|TITLE|RECORD9[A-D]/i and $fields[2] =~ /[^0-9\.+-]/) { warn("The value of $fields[1] ($fields[2]) is not a number - $_\n"); } # now put stuff in the hash: $varname = $fields[1] if($fields[1]); # if there's a name, use it warn("you have not entered a value for $fields[1]\n") if(length($fields[2])==0); if($fields[1]) { # A name was given if(!$P{$fields[1]}) { # ... and it's a new one $P{$varname} = $fields[2]; # put it in the hash } else { # There is a name already $P{$varname} .= "\t$fields[2]"; # append it to any previous values } } else { # this is either a blank name or a name already used. Append to it. $P{$varname} .= "\t$fields[2]"; } } close SHEET; # copy certain chemical data (names in CPRZMvars) to the %P hash: foreach my $v(@CPRZMvars) { $P{$v} = $C{$v} if exists($C{$v}); # reset the ones in the main input file } # Henry's law const has different name and needs a unit conversion (multiply by # 101325 to convert atm to Pa and divide by RT (8.3143*298.15) gives 40.879 $P{HENRYK} = sprintf("%8.2g",$C{HENRY} * 40.879) if exists($C{HENRY}); #$P{HL}=$C{IR}?600:354; # Removed when this part was updated in the #$P{AFIELD}=$C{IR}?172:10; # portion of the script that writes the inp file } ##################################################### sub ReadChem { %Hydrol = (); # remove any hydrolysis values while() { # Second part: read chemical properties chomp; chomp; # go on after a line with the word "Hydrolysis" in it @fields = split(/\t/); # Case: EXAMS commands -- identified by the word EXAMS in field zero if($fields[0] =~ /EXAMS/i) { # Check for any EXAMS commands $fields[1] =~ s/"//g; # remove any " symbols, EXCEL seems to insert them push(@ex, $fields[1]); # when there is a "*" in the field. next; } # skip any (non-EXAMS) lines w/0 data, but keep zeroes: next unless($fields[2] or $fields[2] == '0'); # Case: Hydrolysis -- "pH" appears in the name column-[1] if($fields[1] =~ /pH\s*/i) { $Hydrol{$'} = $fields[2]; # will be converted to rate when solving for the # EXAMS variables KAH, KNH and KBH } # conversion of half-life to rate done in setEXAMSChem # Case: Application interval - name has the word "interval" in it elsif($fields[1] =~ /interval/i and $fields[2] > 0) { # put intervals in their own array push(@AppInt, $fields[2]); # interval is the only multi-valued } # parameter in the chemical file # Otherwise, assume the line is a normal name, value pair: else { $fields[1] =~ s/([a-z])/\U\1/g; # convert name to uppercase for storage $C{$fields[1]} = $fields[2]; # case has already been checked } } delete $Hydrol{""}; # delete any blank key $nApps = @AppInt+1; # this line is because older versions used different values for %C{IR} # so invalid names can be a problem. $C{IR} = "Default" unless($C{IR} =~ /IR|Pond/); # make sure this is valid } ############################################################# sub SetChemParameters { ($C{KAH}, $C{KNH}, $C{KBH}) = CalcHydrolysis(\%Hydrol); # SetEXAMScommands(); # set PRZM degradation rates: # ASM, because it was converted from a half-life, will have more than 8 # decimal places and will mess up the %8g format used to print stuff out my $degrRate = $C{ASM} ? substr(sprintf("%8f", -log(0.5)/$C{ASM}),0,8) : 0; # limit number length $P{DWRATE} = $degrRate; $P{DSRATE} = $degrRate; # these values are the same according to guidance $P{DGRATE} = 0; $P{KD} = $C{KD}; # expand the numbers needed for each horizon into tab-delimited lists: for(my $i=1 ; $i<$P{NHORIZ} ; $i++) { # expand for all horizons $P{DWRATE} .= "\t$degrRate;"; # assumes NHORIZ is known already $P{DSRATE} .= "\t$degrRate;"; $P{DGRATE} .= "\t0"; $P{KD} .= "\t$C{KD}"; } InitPRZMChem(); SetAppSched($C{DATE},@AppInt); } ########################################################### # Sets default values for records 12-18 sub InitPRZMChem { $P{NCHEM} = 1; # for now $P{FRMFLG} = 0; $P{DKFLG2} = 0; } ########################################################### # sets the day and month for each application. It starts with the # date read from the "Date" field of the chemical file and then # adds each interval read from the file to that date. The year # is filled in by the sub PrintR1517 which prints the relevant # records to the PRZM .inp file. sub SetAppSched { my ($date, @AppInt) = @_; my @daysInMonth = (31,28,31,30,31,30,31,31,30,31,30,31); my %months = ("Jan"=>1,"Feb"=>2,"Mar"=>3,"Apr"=>4,"May"=>5,"Jun"=>6, "Jul"=>7,"Aug"=>8,"Sep"=>9,"Oct"=>10,"Nov"=>11,"Dec"=>12); my $intv; my ($d, $m) = split(/\/|-/, $date); if($m =~ /Jan|Feb|Mar|Apr|May|Jun|Jul|Aug|Sep|Oct|Nov|Dec/i) { $m = $months{$m}; } if($m =~ /[^0-9]/ or $m<1 or $m>12) { push @errors, "invalid application month\n"; } ($P{APD}, $P{APM}) = ($d, $m); for my $name("DRFT", "APPEFF","DEPI","TAPP") { $P{$name} = $C{$name}; } foreach $intv(@AppInt) { last if($intv <= 0);# intervals must be greater than zero days $d += $intv; while($d > $daysInMonth[$m-1]) { $d -= $daysInMonth[$m-1]; $m++; push @errors,"Application schedule extends beyond the end of the year\n" if($m>12); } $P{APD} .= "\t$d"; $P{APM} .= "\t$m"; for my $name("DRFT", "APPEFF","DEPI","TAPP") { $P{$name} .= "\t$C{$name}"; } } } ########################################################### sub CalcHydrolysis { my $H = shift @_; # a reference to a hash my $i = 0; my (@hyd, @pH); my ($exKAH, $exKBH, $exKNH); # put hydrolysis information into two parallel arrays and in increasing foreach my $pH(sort keys %$H) { # order of pH. Also substitute 0 half-life $hyd[$i] = $$H{$pH} == 0 ? 0 : log(2)/(24*$$H{$pH}); # for 0 rate $pH[$i++] = $pH; # and convert half-lives (days) into rates (hrs) unless(3 < $pH or $pH < 11) { push @errors, "Is a pH of $pH going to be found in the environment?\n"; } } if(@hyd == 1) { # only one hydrolysis value entered use it for KNH $exKNH = $hyd[0]; $exKBH = 0; $exKAH = 0; } elsif(@hyd == 3) { ($exKAH, $exKNH, $exKBH) = HydSolve(scalar(@hyd), \@pH, \@hyd); } else { push @errors, "Can't run with ",scalar(@hyd)," hydrolysis values\n"; } return ($exKAH, $exKNH, $exKBH); } ############################################### sub HydSolve { # solve for EXAMS hydrolysis constants my ($n, $pHr, $hydr) = @_; my @A; die("Currently HydSolve only works with three data pairs\n") if($n != 3); for(my $i=0 ; $i<$n ; $i++) { $A[$i][0] = 10.0**(-$$pHr[$i]); $A[$i][1] = 1; $A[$i][2] = 10**($$pHr[$i] - 14); } # PA(\@A, $n); my @x = (0) x $n; my @y = (0) x $n; LUdecompose(\@A, $n); # forward substitution with L for(my $j=0 ; $j < $n ; $j++) { $y[$j] = $$hydr[$j]; for(my $i=0 ; $i<$j ; $i++) { $y[$j] -= $A[$j][$i]*$y[$i]; } # Would divide by A[j][j], but that's 1 because we're using L # print "y[$j] = $y[$j]\n"; } # back substitution with U to obtain the answer for(my $j=$n-1 ; $j>=0 ; $j--) { $x[$j] = $y[$j]; for(my $i=$n-1 ; $i>$j ; $i--) { $x[$j] -= $A[$j][$i]*$x[$i]; } $x[$j] /= $A[$j][$j]; # print "x[$j] = ",$x[$j],"\n"; } return @x; } ############################################################# # Here we do gaussian elimination to decompose the matrix. The part # below the diagonal is L-I (I=identity matrix) and the diagonal and # above is U. sub LUdecompose { my ($a, $n) = @_; # a reference to the matrix and the dimension for my $i(0..$n-2) { for my $j($i+1..$n-1) { $a->[$j][$i] /= $a->[$i][$i]; } for my $j($i+1..$n-1) { for my $k($i+1..$n-1) { $a->[$j][$k] -= $a->[$j][$i] * $a->[$i][$k]; } } } } ############################################################ sub PA { my ($A, $n) = @_; for($i=0 ; $i<$n ; $i++) { for($j=0 ; $j<$n ; $j++) { printf "%10g\t", $A->[$i][$j]; } print "\n"; } print "\n"; } ######################## # This does just what it says: it extracts the dates from the # first and last lines of the metfile # It could also write a metfile from an EXPRES .AWD file sub findmetbounds { debug("+findmetbounds -- $metfilename"); if($metfilename !~ /\.dvf|\.met|\.awd/i) { push @errors, ("metfiles should have the extension .dvf, .met or .awd\n"); } my ($first, $last); my ($mnth, $wind, $temp); my @days = (); @Tmp_mnth = (); @Wnd_mnth = (); if($metfilename =~ /\.awd/i) { ($metfilename,$first,$last) = ConvertEXPRESmetfile($`); # send the name without the extension } open(MET, $filepath."\\metfiles\\".$metfilename) or push @errors, "can't open metfile $metfilename\n"; $_ = ; $first = substr($_,1,6); $mnth = substr($_,1,2) - 1; # zero offset months go from 0 to 11 my $yr1 = substr($_,5,2); my $yr = 0; $days[0][$mnth]++; $Tmp_mnth[0][$mnth] += substr($_,27,10); $Wnd_mnth[0][$mnth] += substr($_,37,10); $last = "error"; while(){ $last = substr($_,1,6); $mnth = substr($_,1,2) - 1; # zero offset months go from 0 to 11 $yr = substr($_,5,2) - $yr1; $days[$yr][$mnth]++; $Tmp_mnth[$yr][$mnth] += substr($_,27,10); $Wnd_mnth[$yr][$mnth] += substr($_,37,10); } close MET; push @errors, "error reading metfile\n" if($last eq "error"); for($yr = 0 ; $yr < @days ; $yr++) { for $mnth(0..$#{$days[$yr]}) { my $d = $days[$yr][$mnth]; if($d > 0) { $Tmp_mnth[$yr][$mnth] /= $d; $Wnd_mnth[$yr][$mnth] /= $d; } else { $Tmp_mnth[$yr][$mnth] = "No Data"; $Wnd_mnth[$yr][$mnth] = "No Data"; } } } return ($first, $last); } ############################################ # Converts an .awd file from EXPRES into PRZM metfile format. This routine # generates a new file with the .met extension (and the original filename) # and leaves it in the metfiles directory. sub ConvertEXPRESmetfile { my $name = shift; $namepath = $filepath."\\metfiles\\".$name; print "metfile with path is $namepath\n"; open AWD, "$namepath.awd" or die("can't open EXPRES weather file $name.awd \n"); open MET, ">$namepath.met" or die("can't open MET file $name.met - $?\n"); $_ = ; # skip one $_ = ; # get initial date from the second line chomp; chomp; my $f = substr($_,2,2).substr($_,0,2).substr($_,4,2); my $l = "error"; my $z = 0; printf MET " %02d%02d%02d%10g%10g%10g%10g%10g\n",substr($_,2,2),substr($_,0,2), substr($_,4,2),substr($_,6,5)/10,substr($_,11,6)/10, substr($_,17,6), 0,0; # the AWD file has no data for windspeed or solar radiation while() { #cycle through the remaining lines chomp; chomp; last if(length() < 20); #This line is too short to contain met data, it must be the end. printf MET " %02d%02d%02d%10g%10g%10g%10g%10g\n",substr($_,2,2),substr($_,0,2), substr($_,4,2),substr($_,6,5)/10,substr($_,11,6)/10, substr($_,17,6), 0,0; # the AWD file has no data for windspeed or solar radiation $l = substr($_,0,6); } close AWD; close MET; $l = substr($l,2,2).substr($l,0,2).substr($l,4,2); return ("$name.met", $f, $l); } ############################################ sub TestData { my $compSum = 0; my $nCompart; my $d = 0; my $i = 0; my @dpn = split(/\t/, $P{DPN}); my @warnings = (); foreach $d(split(/\t/, $P{THKNS})) { my $test = $d / $dpn[$i]; if(($test - int($test)) > 0.0001) { push @warnings,"non-integer number of compartments in layer $i"; } $compSum += $d; # sum up the total thickness of the layers $nCompart += $d/$dpn[$i++]; # sum up the number of compartments } if($compSum != ($P{CORED})) { push @warnings, "The sum of the soil layers ($compSum) != core depth" . "($P{CORED}) -- resetting CORED to $compSum."; $P{CORED} = $compartmentsum; } push @warnings, "There are more than 500 ($nCompart) ". "compartments, exceeding PRZM's array size" if($nCompart > 500); my $rootdep = 0; foreach $d (split(/\t/, $P{AMXDR})) { $rootdep = $d if($d > $rootdep); } if($rootdep > $P{CORED}) { push @warnings, "maximum rooting depth ($rootdep) below bottom of soil profile $P{CORED}"; } $P{NDC} = 1 if($P{NDC} < 1); push @warnings, "You have $P{NDC} crops\n" if($P{NDC} > 1); my $result = "Continue"; # Check if "EXAMS.daf" is available: No need for EXAMS 2.98 # if(opendir DIR, "./") { # my @files = grep /EXAMS\.daf/i, readdir(DIR); # push @warnings, "EXAMS.DAF cannot be found." unless @files; # } # else { # push @warnings, "Can't open current directory."; # } if($EXAMSenv =~ /IR|Res/i and $C{IR} !~ /IR/i) { push @warnings, "You are running EXAMS with an Index Reservoir. Check PRZM field area."; } elsif($EXAMSenv =~ /Pond/i and $C{IR} !~ /pond/i) { push @warnings, "You are running EXAMS with the pond. Check PRZM field area."; } if($EXAMSenv =~ /IR|Res/i and $C{RUNOFF} =~ /none/i) { push @warnings, "You might want to use runoff"; } elsif($EXAMSenv =~ /pond/i and $C{RUNOFF}=~ /total|monthly/i) { push @warnings, "Are you sure you want flow through the pond?"; } if(@warnings) { my $warnDial = $mainwindow->DialogBox(-title => "PRZM parameter warnings", -buttons => ["Continue", "Abort", "Exit"]); my $warnBox = $warnDial->Listbox( -width => 60, -height => 10)-> pack(-side => "top", -anchor => "nw"); foreach (@warnings) { $warnBox->insert('end', $_); } my $scroll = $warnDial->Scrollbar( -command => ['yview', $warnBox]); $warnBox->configure(-yscrollcommand => ['set', $scroll]); $warnBox->pack(-side => 'left', -fill => 'both', -expand => 'yes'); $scroll->pack(-side => 'right', -fill => 'y'); $result = $warnDial->Show; # Execute the dialog box } exit if $result eq "Exit"; return $result eq "Continue" ? 1 : 0; } ################################################ # Estimates the DAIR parameter based on the diffusion coefficient of diethyl # ether and the sqrt of the ratio of molecular weights. sub DAIRest { my $molWt = shift; if($C{MWT} <= 0) { # this is an error: no molecular weight # set some error message return 0; } my $D = 74/$molWt; $D = sqrt($D); $D *= 7689.6; # diff. coeff of diethyl ether in cm^2/day $D = sprintf("%7f",$D); return $D; } ############################################ # Generate the PRZM input file: sub WritePRZMinput { debug("+WritePRZMinput"); # first set any values of the %P hash stored in other variables: foreach my $v(@CPRZMvars) { $P{$v} = $C{$v} if exists($C{$v}); # reset the ones in the main input file } # Set a few variables which we don't set explicitly earlier. These are based on # the values of other variables already entered. # Henry's law const has different name and needs a unit conversion (multiply by # 101325 to convert atm to Pa and divide by RT (8.3143*298.15) gives 40.879 $P{HENRYK} = sprintf("%8.2g",$C{HENRY} * 40.879) if exists($C{HENRY}); $P{PTITLE} = "$P{PSTNAM} - $apps applications @ $P{TAPP} kg/ha"; if($C{KD} <= 0) { if($C{KOC} > 0) { $P{KDFLAG} = 1; $P{PCMC} = 4; $P{SOL} = $C{KOC}; } elsif($C{SOL} > 0) { $P{KDFLAG} = 1; $P{PCMC} = 2; # mg/L $P{SOL} = $C{SOL}; } } else { $P{KDFLAG} = 0; } # Start writing the input file open(PINP, ">przm3.inp") || die("Can't open przm3.inp for writing\n"); printf PINP "%-78s\n", $P{TITLE}; printf PINP "%-78s\n", $P{HTITLE}; print PINP "*** Record 3:\n"; printf PINP "%8g%8g%8d%8g%8d%8d\n", $P{PFAC}, $P{SFAC}, $P{IPEIND}, $P{ANETD}, $P{INICRP}, $P{ISCOND}; # we're not using WDM data, so DSN is not printed out if($P{IPEIND} == 1 || $P{IPEIND} == 2 ) { # Records 4 and 5 print PINP "*** Record 4 (followed by Record 5)\n"; my @DT = split(/\t/, $P{DT}); for(my $i=0 ; $i<12 ; $i++) { printf PINP "%8g", $DT[$i]; print(PINP "\n") if($i%6 == 5); # put a new line after every 6 values for DT } } print PINP "*** Record 6 -- ERFLAG\n"; printf PINP "%8d\n", $P{ERFLAG}; if(2 <= $P{ERFLAG} && $P{ERFLAG} <= 4) { print PINP "*** Record 7:\n"; my ($afield, $hl) = FieldSize(); printf PINP "%8g%8g%8g%8g %8d%8g%8g\n", $P{USLEK}, $P{USLELS}, $P{USLEP}, $afield, $P{IREG}, $P{SLP}, $hl; } print PINP "*** Record 8\n"; printf PINP "%8d\n", $P{NDC}; PrintRecord9(); # prints 9, 9A, 9B, 9C and 9D for all crops. print PINP "*** Record 10 -- NCPDS, the number of cropping periods\n"; printf PINP "%8d\n", $lastYr - $firstYr + 1; PrintRecord11(); print PINP "*** Record 12 -- PTITLE\n"; printf PINP "%-78s\n", $P{PTITLE}; print PINP "*** Record 13\n"; printf PINP "%8d%8d%8d%8d\n", $P{NAPS}, $P{NCHEM}, $P{FRMFLG}, $P{DKFLG2}; if($P{DKFLG2} == 1) { my @DKDAY = split(/\t/, $P{DKDAY}); my @DKMNTH = split(/\t/, $P{DKMNTH}); my @DKNUM = split(/\t/, $P{DKNUM}); if(scalar(@DKDAY) != $P{NCHEM} || scalar(@DKMNTH) != $P{NCHEM} || scalar(@DKNUM) != $P{DCNUM}) { die("Number of parameters does not match NCHEM in R14\n") ; } for(my $i=0 ; $i<$P{NCHEM} ; $i++) { printf PINP " %2.2d%2.2d%8d", $DKDAY[$i], $DKMNTH[$i], $DKNUM[$i]; } print PINP "\n"; } # Return value of Print$1517,$cam23, tells if CAM ever has a value of 2 or 3 if(PrintR1517()) { print PINP "*** Record 18\n"; my @PLVKRT = split(/\t/, $P{PLVKRT}); my @PLDKRT = split(/\t/, $P{PLDKRT}); my @FEXTRC = split(/\t/, $P{FEXTRC}); for(my $i=0 ; $i<$P{NCHEM} ; $i++) { printf PINP "%8g%8g%8g", $PLVKRT[$i], $PLDKRT[$i], $FEXTRC[$i]; } print PINP "\n"; if($P{NCHEM} > 1) { print PINP "*** Record 18A\n"; my @PTRAN12 = split(/\t/, $P{PTRAN12}); my @PTRAN13 = split(/\t/, $P{PTRAN13}); my @PTRAN23 = split(/\t/, $P{PTRAN23}); printf PINP "%8g%8g%8g\n", $P{PTRAN12}, $P{PTRAN13}, $P{PTRAN23}; } } print PINP "*** Record 19 -- STITLE\n"; printf PINP "%-78s\n", $P{STITLE}; print PINP "*** Record 20\n"; printf PINP "%8g %4d%4d%4d%4d%4d%4d%4d%4d%4d\n", $P{CORED}, $P{BDFLAG}, $P{THFLAG}, $P{KDFLAG}, $P{HSWZT}, $P{MOC}, $P{IRFLAG}, $P{ITFLAG}, $P{IDFLAG}, $P{BIOFLG}; if($P{BIOFLG} == 1) { print PINP "*** Record 21\n"; printf PINP "%8g%8g%8g%8g%8g\n", $P{AM}, $P{AC}, $P{AS}, $P{AR}, $P{KE}; print PINP "*** Record 22\n"; printf PINP "%8g%8g%8g%8g%8g%8g%8g\n", $P{KSM}, $P{KCM}, $P{KC}, $P{MKS}, $P{KR}, $P{KIN}, $P{KSK}; print PINP "*** Record 23\n"; printf PINP "%8g%8g%8g%8g%8g%8g%\n", $P{KLDM}, $P{KLDC}, $P{KLDS}, $P{KLDR}, $P{KL1}, $P{KL2}; print PINP "*** Record 24\n"; printf PINP "%8g%8g%8g%8g%8g\n", $P{USM}, $P{UCM}, $P{MUC}, $P{US}, $P{UR}; print PINP "*** Record 25\n"; printf PINP "%8g%8g%8g%8g%8g\n", $P{YSM}, $P{YCM}, $P{YC}, $P{US}, $P{YR}; } print PINP "*** Record 26\n"; { # This block is just to make the "my" variables below expire when they're not needed my @dair = split(/\t/, $P{DAIR}); my @henryk = split(/\t/, $P{HENRYK}); my @enpy = split(/\t/, $P{ENPY}); for(my $i=0 ; $i<$P{NCHEM} ; $i++) { unless($enpy[$i]) { $dair[$i] = 0; $henryk[$i] = 0; } elsif($henryk[$i] != 0 and $dair[$i] == 0) { $dair[$i] = DAIRest($C{MWT}); # We would need more than one MWT for multiple chems } # note printf seems to use 3-digit exponents, so the decimals of HENRYK are limited printf PINP "%8g%8.2g%8g", $dair[$i], $henryk[$i], $enpy[$i]; } print PINP "\n"; } # Records 27-29 - Irrigation parameters if($P{IRFLAG} == 1 or $P{IRFLAG} == 2) { print PINP "*** Record 27 -- irrigation\n"; printf PINP "%8d%8g%8g%8g\n", $P{IRTYP}, $P{FLEACH}, $P{PCDEPL}, $P{RATEAP}; if($P{IRTYP} == 2) { print PINP "*** Record 28 - Furrow irrigation paramters\n"; printf PINP "%8g%8g%8g%8g%8g%8g%8g\n", $P{Q0}, $P{BT}, $P{ZRS}, $P{SF}, $P{EN}, $P{X2}, $P{XFRAC}; print PINP "*** Record 29 - more furrow irrigation parameters\n"; printf PINP "%8g%8g\n", $P{KS}, $P{HF}; } } if($P{KDFLAG} == 1) { print PINP "*** Record 30\n"; printf PINP "%8d", $P{PCMC}; my @sol = split(/\t/, $P{SOL}); for(my $i=0 ; $i<$P{NCHEM} ; $i++) { printf PINP "%8g", $sol[$i]; } print PINP "\n"; } if($P{ITFLAG} == 1) { print PINP "*** Record 31\n"; my @ALBEDO = split(/\t/, $P{ALBEDO}); die("Wrong number of albedo values, R31\n") if(scalar(@ALBEDO) != 12); for(my $i = 0 ; $i<12 ; $i++) { printf PINP "%5g", $ALBEDO[$i]; } printf PINP "%5g%5g\n", $P{EMMISS}, $P{ZWIND}; my @BBT = split(/\t/, $P{BBT}); die("Wrong number of albedo values, R31\n") if(scalar(@BBT) != 12); for(my $i = 0 ; $i<12 ; $i++) { printf PINP "%5g", $BBT[$i]; } printf PINP "%8g%8g\n", $P{QFAC}, $P{TBASE}; } PrintSoil(); printf PINP "***Record 40\n%8d", $P{ILP}; if($P{ILP} == 1) { printf PINP "%8d\n", $P{CLFAG}; print "Error: this script does not currently print out record 41.\n"; die(" ILP should be set to zero\n"); } else {print PINP "\n"}; # Print records 42 and higher. This part is just inserted verbatim from mscot.inp. print PINP " YEAR 10 YEAR 10 YEAR 10 1\n"; print PINP " 1\n"; print PINP " 1 -----\n"; print PINP " 7 YEAR\n"; print PINP " PRCP TCUM 0 0\n"; print PINP " RUNF TCUM 0 0\n"; print PINP " INFL TCUM 1 1\n"; print PINP " ESLS TCUM 0 0 1.0E3\n"; print PINP " RFLX TCUM 0 0 1.0E5\n"; print PINP " EFLX TCUM 0 0 1.0E5\n"; print PINP " RZFX TCUM 0 0 1.0E5\n"; close PINP; debug("-WritePRZMinput"); } ############ END OF PRZM INPUT FILE ############### ############################# sub PrintRecord9 { print PINP "*** Record 9\n"; my @ICNCN = split(/\t/, $P{ICNCN}); my @CINTCP = split(/\t/, $P{CINTCP}); my @AMXDR = split(/\t/, $P{AMXDR}); my @COVMAX = split(/\t/, $P{COVMAX}); my @ICNAH = split(/\t/, $P{ICNAH}); my @CN = split(/\t/, $P{CN}); my @WFMAX; if($P{CAM} =~ /3/) { @WFMAX = split(/\t/, $P{WFMAX}); } else { @WFMAX = (0) x scalar(@ICNAH); } my @HTMAX = split(/\t/, $P{HTMAX}); for(my $i = 0 ; $i<$P{NDC} ; $i++) { printf PINP "%8d%8g%8g%8g%8d %3d %3d %3d%8g%8g\n", $ICNCN[$i], $CINTCP[$i], $AMXDR[$i], $COVMAX[$i], $ICNAH[$i], $CN[3*$i], $CN[3*$i+1], $CN[3*$i+2], $WFMAX[$i], $HTMAX[$i]; } if(exists $P{RECORD9A} and exists $P{RECORD9B} and exists $P{RECORD9C} and exists $P{RECORD9D}) { print PINP "*** Record 9a-d\n"; my @R9A = split(/\t/, $P{RECORD9A}); my @R9B = split(/\t/, $P{RECORD9B}); my @R9C = split(/\t/, $P{RECORD9C}); my @R9D = split(/\t/, $P{RECORD9D}); my $ndc = scalar(@R9A); my $nlines = scalar(@R9B)/$ndc; for my $j(0..$ndc-1) { print PINP "$R9A[$j]\n"; for my $i(0..$nlines-1) { print PINP $R9B[$j*2+$i], "\n"; print PINP $R9C[$j*2+$i], "\n"; print PINP $R9D[$j*2+$i], "\n"; } } } else { my @CROPNO = split(/\t/, $P{CROPNO}); my @NUSLEC = split(/\t/, $P{NUSLEC}); my @GDUSLEC = split(/\t/, $P{GDUSLEC}); my @GMUSLEC = split(/\t/, $P{GMUSLEC}); my @USLEC = split(/\t/, $P{USLEC}); my @MNGN = split(/\t/, $P{MNGN}); my $uslec1 = 0; # the index of the first USLEC and MNGN of the current crop # it is incremented as we cycle through each crop. for(my $i = 0 ; $i<$P{NDC} ; $i++) { printf PINP "%8d%8d\n", $CROPNO[$i], $NUSLEC[$i]; # 9A my $mx = ($NUSLEC[$i] <= 16) ? $NUSLEC[$i] : 16; for(my $j = 0 ; $j<$mx; $j++) { # 9B printf PINP "%2.2d%2.2d ", $GDUSLEC[$uslec1+$j], $GMUSLEC[$uslec1+$j]; } print PINP "\n"; for(my $j = 0 ; $j<$mx ; $j++) { # 9C printf PINP "%4.2f ", $USLEC[$uslec1+$j]; } print PINP "\n"; for(my $j = 0 ; $j<$mx ; $j++) { # 9D printf PINP "%4g ", $MNGN[$uslec1+$j]; } print PINP "\n"; # The following section prints Records 9BCD if there are # more than 16 USLEC and MNGN numbers. if($NUSLEC[$i] > 16 ) { for(my $j = $mx ; $j<$NUSLEC[$i] ; $j++) { # 9B printf PINP "%2.2d%2.2d ", $GDUSLEC[$uslec1+$j], $GMUSLEC[$uslec1+$j]; } print PINP "\n"; for(my $j = $mx ; $j<$NUSLEC[$i]; $j++) { # 9C printf PINP "%4g ", $USLEC[$uslec1+$j]; } print PINP "\n"; for(my $j = $mx ; $j<$NUSLEC[$i]; $j++) { # 9D printf PINP "%4g ", $MNGN[$uslec1+$j]; } print PINP "\n"; } $uslec1 += $NUSLEC[$i]; } } # end if } ################################### # Subroutine for printing Record 11, crop emergence, maturity and harvest # times. # Right now this routine ignores the year values (IYREM, IYRMAT, and IYRHAR) # and instead uses the years in the metfile, cycling through the list of # crops. sub PrintRecord11 { print PINP "*** Record 11\n"; my @emd = split(/\t/, $P{EMD}); my @emm = split(/\t/, $P{EMM}); my @mad = split(/\t/, $P{MAD}); my @mam = split(/\t/, $P{MAM}); my @had = split(/\t/, $P{HAD}); my @ham = split(/\t/, $P{HAM}); my @emy = split(/\t/, $P{IYREM}); my @may = split(/\t/, $P{IYRMAT}); my @hay = split(/\t/, $P{IYRHAR}); my (@yrOffsetMat, @yrOffsetEm); for(my $i=0 ; $i<@emy ; $i++) { $yrOffsetMat[$i] = $may[$i] - $hay[$i]; $yrOffsetEm[$i] = $emy[$i] - $hay[$i]; } my @incrop = split(/\t/, $P{INCROP}); for(my $i=$firstYr ; $i<=$lastYr ; $i+=$P{NDC}) { for(my $j=0 ; $j<$P{NDC} ; $j++) { last if($i+$j > $lastYr); printf PINP " %2.2d%2.2d%2.2d %2.2d%2.2d%2.2d %2.2d%2.2d%2.2d%8d\n", $emd[$j], $emm[$j], $j + $i + $yrOffsetEm[$j], $mad[$j], $mam[$j], $j + $i + $yrOffsetMat[$j], $had[$j], $ham[$j], $j + $i, $incrop[$j]; } } } #################################### # prints Pesticide names and application schedule. # Currently this only works with one application date/year. sub PrintR1517 { print PINP "*** Record 15 -- PSTNAM\n"; # split doesn't seem to work if there is nothing to split on # for now just spit out the PSTNAM parameter on this line # my $PSTNAM = split(/\t/, $P{PSTNAM}); # for($i=0 ; $i<$P{NCHEM} ; $i++) { # printf PINP "%-20s", $PSTNAM[$i]; # } print PINP $P{PSTNAM}; print PINP "\n*** Record 16\n"; my @apd = split(/\t/, $P{APD}); my @apm = split(/\t/, $P{APM}); my @winday = split(/\t/, $P{WINDAY}); my @cam = split(/\t/, $P{CAM}); my @depi = split(/\t/, $P{DEPI}); my @tapp = split(/\t/, $P{TAPP}); my @appeff = split(/\t/, $P{APPEFF}); my @drft = split(/\t/, $P{DRFT}); my $cam23 = 0; # flag to tell if any value of CAM is equal to 2 or 3 for(my $y = $firstYr ; $y<=$lastYr ; $y++) { # yearly loop for(my $i=0 ; $i<$nApps ; $i++) { #loop for apps within each year printf PINP " %2.2d%2.2d%2.2d%3d", $apd[$i], $apm[$i], $y, $winday[$i]; for(my $j=0 ; $j<$P{NCHEM} ; $j++) { printf PINP "%2d%5.1f%6g%5.3g%5.3g", $cam[$j], $depi[$j], $tapp[$j], $appeff[$j], $drft[$j]; $cam23 = 1 if($cam[$j] == 2 or $cam[$j] == 3); } print PINP "\n"; } } my @IPSCND = split(/\t/, $P{IPSCND}); my @UPTKF = split(/\t/, $P{UPTKF}); print PINP "*** Record 17\n"; $P{FILTRA} = 0 unless($cam[0] == 3); printf PINP "%8g", $P{FILTRA}; for(my $i=0 ; $i<$P{NCHEM} ; $i++ ) { printf PINP "%8d%8g", $IPSCND[$i], $UPTKF[$i]; } print PINP "\n"; return $cam23; } ##################################### sub PrintSoil { # NOTE: This routine does not currently set record 39 for cases when DKFLG2 # is set (bi-phasic half-life). print PINP "*** Record 33\n"; printf PINP "%8d\n", $P{NHORIZ}; my @HORIZN = split(/\t/, $P{HORIZN}); my @THKNS = split(/\t/, $P{THKNS}); my @BD = split(/\t/, $P{BD}); my @THETO = split(/\t/, $P{THETO}); my @AD = split(/\t/, $P{AD}); my @DISP = split(/\t/, $P{DISP}); my @ADL = split(/\t/, $P{ADL}); my @DPN = split(/\t/, $P{DPN}); my @THEFC = split(/\t/, $P{THEFC}); my @THEWP = split(/\t/, $P{THEWP}); my @OC = split(/\t/, $P{OC}); my @KD = split(/\t/, $P{KD}); my @AD = split(/\t/, $P{AD}); my $n = $P{NCHEM}; # just to make it shorter for(my $i=0 ; $i<$P{NHORIZ} ; $i++) { # print PINP "*** Record 34\n"; printf PINP "%8d%8g%8g%8g%8g%8g%8g\n", $HORIZN[$i], $THKNS[$i], $BD[$i], $THETO[$i], $AD[$i], $DISP[$i], $ADL[$i]; if($P{BIOFLG} == 1) { # Record 35 my @Q = split(/\t/, $P{Q}); my @CM1 = split(/\t/, $P{CM1}); my @Y1 = split(/\t/, $P{Y1}); my @Y2 = split(/\t/, $P{Y2}); my @Y3 = split(/\t/, $P{Y3}); my @Y4 = split(/\t/, $P{Y4}); printf PINP " %8g%8g%8g%8g%8g\n", $Q[$i], $CM1[$i], $Y1[$i], $Y2[$i], $Y3[$i], $Y4[$i]; } if($P{DKFLG2} == 0) { # Record 36 my @DWRATE = split(/\t/, $P{DWRATE}); my @DSRATE = split(/\t/, $P{DSRATE}); my @DGRATE = split(/\t/, $P{DGRATE}); print PINP " "; for(my $j=0 ; $j< $P{NCHEM} ; $j++ ) { printf PINP "%8g%8g%8g", $DWRATE[$i*$n+$j], $DSRATE[$i*$n+$j], $DGRATE[$i*$n+$j]; } print PINP "\n"; } if($P{DKFLG2} == 1) { my @DWRAT1 = split(/\t/, $P{DWRAT1}); my @DSRAT1 = split(/\t/, $P{DSRAT1}); my @DGRAT1 = split(/\t/, $P{DGRAT1}); my @DWRAT2 = split(/\t/, $P{DWRAT2}); my @DSRAT2 = split(/\t/, $P{DSRAT2}); my @DGRAT2 = split(/\t/, $P{DGRAT2}); print PINP " "; for(my $j=0 ; $j< $n ; $j++ ) { printf PINP "%8g%8g%8g", $DWRAT1[$i*$n+$j], $DSRAT1[$i*$n+$j], $DGRAT1[$i*$n+$j]; } print PINP "/n"; print PINP " "; for(my $j=0 ; $j< $n ; $j++ ) { printf PINP "%8g%8g%8g", $DWRAT2[$i*$n+$j], $DSRAT2[$i*$n+$j], $DGRAT2[$i*$n+$j]; } print PINP "/n"; } # print PINP "*** Record 37\n"; printf PINP " %8g%8g%8g%8g", $DPN[$i], $THEFC[$i], $THEWP[$i], $OC[$i]; for(my $j=0 ; $j< $P{NCHEM} ; $j++ ) { printf PINP "%8g", $KD[$i*$n+$j]; } print PINP "\n"; if($P{ITFLAG} == 1) { my @SPT = split(/\t/, $P{SPT}); my @SAND = split(/\t/, $P{SAND}); my @CLAY = split(/\t/, $P{CLAY}); my @THCOND = split(/\t/, $P{THCOND}); my @VHTCAP = split(/\t/, $P{VHTCAP}); printf PINP " %8g%8g%8g", $SPT[$i], $SAND[$i], $CLAY[$i]; if($P{IDFLAG} == 0) {printf PINP "%8g%8g\n", $THCOND[$i], $VHTCAP[$i];} else { printf PINP " 0.0 0.0\n";} } if($P{DKFLG2} == 0 && $n > 1) { print PINP "*** Record 39\n"; my @DKRW12 = split(/\t/, $P{DKRW12}); my @DKRS12 = split(/\t/, $P{DKRS12}); my @DKRW13 = $P{NCHEM} == 2 ? (0.0)x $P{NHORIZ} : split(/\t/, $P{DKRW13}); my @DKRW23 = $P{NCHEM} == 2 ? (0.0)x $P{NHORIZ} : split(/\t/, $P{DKRW23}); my @DKRS13 = $P{NCHEM} == 2 ? (0.0)x $P{NHORIZ} : split(/\t/, $P{DKRS13}); my @DKRS23 = $P{NCHEM} == 2 ? (0.0)x $P{NHORIZ} : split(/\t/, $P{DKRS23}); printf PINP " %8g%8g%8g%8g%8g%8g\n", $DKRW12[$i], $DKRW13[$i], $DKRW23[$i], $DKRS12[$i], $DKRS13[$i], $DKRS23[$i]; } if($P{DKFLG2} == 1 && $n > 1) { my @DKW112 = split(/\t/, $P{DKW112}); my @DKW113 = $P{NCHEM} == 2 ? (0.0)x $P{NHORIZ} : split(/\t/, $P{DKW113}); my @DKW123 = $P{NCHEM} == 2 ? (0.0)x $P{NHORIZ} : split(/\t/, $P{DKW123}); my @DKS112 = split(/\t/, $P{DKS112}); my @DKS113 = $P{NCHEM} == 2 ? (0.0)x $P{NHORIZ} : split(/\t/, $P{DKS113}); my @DKS123 = $P{NCHEM} == 2 ? (0.0)x $P{NHORIZ} : split(/\t/, $P{DKS123}); my @DKW212 = split(/\t/, $P{DKW212}); my @DKW213 = $P{NCHEM} == 2 ? (0.0)x $P{NHORIZ} : split(/\t/, $P{DKW213}); my @DKW223 = $P{NCHEM} == 2 ? (0.0)x $P{NHORIZ} : split(/\t/, $P{DKW223}); my @DKS212 = split(/\t/, $P{DKS212}); my @DKS213 = $P{NCHEM} == 2 ? (0.0)x $P{NHORIZ} : split(/\t/, $P{DKS213}); my @DKS223 = $P{NCHEM} == 2 ? (0.0)x $P{NHORIZ} : split(/\t/, $P{DKS223}); print PINP "*** Record 39\n"; printf PINP " %8g%8g%8g%8g%8g%8g\n", $DKW112[$i], $DKW113[$i], $DKW123[$i], $DKS112[$i], $DKS113[$i], $DKS123[$i]; print PINP "*** Record 39A\n"; printf PINP " %8g%8g%8g%8g%8g%8g\n", $DKW212[$i], $DKW213[$i], $DKW223[$i], $DKS212[$i], $DKS213[$i], $DKS223[$i]; } } return; } # End of SoilPrint; ################################################################### # returns values for $afield and $hl based on the setting of $C{IR} sub FieldSize { my ($afield, $hl) = ($P{AFIELD}, $P{HL}); if($C{IR} eq "IR"){ ($afield, $hl) = (172.8, 600); } elsif($C{IR} eq "Pond") { ($afield, $hl) = (10, 354); } return ($afield, $hl); } ################################################################### sub WriteRunFile { open(PRUN, ">przm3.run") or die "couldn't open przm3.run for writing\n"; print PRUN "*** Options\n"; print PRUN "PRZM ON\n"; print PRUN "VADOFT OFF\n"; print PRUN "MONTE CARLO OFF\n"; print PRUN "TRANSPORT OFF\n"; print PRUN "*** Zone records\n"; print PRUN "PRZM ZONES 1\n"; print PRUN "ENDRUN\n"; print PRUN "*** input file records\n"; # print PRUN " PATH $path\\\n"; # metfile has a different path print PRUN " METEOROLOGY 1 ",$filepath."\\metfiles\\".$metfilename,"\n"; warn "Generating runfile: metfilename is $metfilename\n"; print PRUN " PRZM INPUT 1 przm3.inp\n"; print PRUN "*** output file records\n"; print PRUN " PATH $path\\\n"; # ($file) = split(/\./, $file); # shave the extension off the input spreadsheet name print PRUN " TIME SERIES 1 $outputname".".zts\n"; print PRUN " PRZM OUTPUT 1 $outputname".".zpm\n"; print PRUN "*** scratch file records\n"; print PRUN " PATH $path\\\n"; print PRUN " PRZM RESTART RESTART.PRZ\n"; print PRUN "ENDFILES\n"; print PRUN "*** global records\n"; my $sd = substr($startdate,2,2).substr($startdate,0,2).substr($startdate,4,2); print PRUN " START DATE $sd\n"; my $ed = substr($enddate,2,2).substr($enddate,0,2).substr($enddate,4,2); print PRUN " END DATE $ed\n"; printf PRUN "NUMBER OF CHEMICALS %-2d\n", $P{NCHEM}; print PRUN "ENDDATA\n"; print PRUN "*** display records\n"; print PRUN "ECHO 8\n"; print PRUN "TRACE ON\n"; close PRUN; } ########################################################### # Take the values from the chemical sheet and write out commands to EXAMS # to set them in that model. This stores those commands in the array @ex. # They are written to the .exa file later in the sub WriteEXAfile. # This uses unshift to add commands to the array to be sure the chemical settings # come first in the .exa file, user settings, added with push, will be at the # end, no matter which order they were added in. sub SetEXAMSChem { print EXA "recall chem 1\n"; print EXA "chemical name is $P{PSTNAM}\n"; # KD is used by both, %C has it as single value print EXA "set KPS(1,1) to $C{KD}\n"if($C{KD} ne ""); my ($s, $var, $name, $value); my @toset = ("MWT(1)", "HENRY(1)", "VAPR(1)", "SOL(1,1)", "KDP(1,1)", "KBACW(*,1,1)", "KBACS(*,1,1)", "KAH(1,1,1)", "KNH(1,1,1)", "KBH(1,1,1)", "KOC(1)"); foreach $var(@toset) { $var =~ /\(/; # this is to get the name part of var isolated in $` if($C{$`}) { if($var =~ /(KDP|KBAC.)/ ) { # convert these to rates $value = $C{$1} == 0 ? 0 : log(2)/(24*$C{$1}); # t1/2 of 0 means rate of 0 } else { $value = $C{$`}; } print EXA "set $var = $value\n" if($value ne ""); } } # the Q10 parameters are set constant in this script: print EXA "set QTBAS(*,1,1) = 2\n"; print EXA "set QTBAW(*,1,1) = 2\n"; } ############################################### # Writes the EXAMS batch file sub WriteEXAfile { opendir(DIR, "./") or die("Couldn't open current directory ./\n"); @files = grep /P2E-/, readdir DIR; open(EXA, ">pz2ex.exa") or die("Can't open pz2ex.exa for writing\n"); print EXA "set mode = 3\n"; print EXA "set outfil(4) to Y\n"; print EXA "set outfil(2) to N\n"; print EXA "READ ENV ",$filepath."\\EXAMSenv\\".$EXAMSenv,"\n"; print EXA "READ MET ",$filepath."\\Metfiles\\".$metfilename,"\n"; print EXA "SET YEAR1 = ",$firstYr < 100 ? $firstYr+1900 : $firstYr,"\n"; SetEXAMSChem(); foreach $cmd(@ex) { print EXA "$cmd\n"; } my $rf; # stores total average runoff for the run, if needed if($C{RUNOFF} =~ /total/i) { # monthly runoff is handled in the loop below $rf = GetTotalRunoff(); } my $run = 0; # a flag my $yr = 0; # $offset sets how many months before the current month to use temps from my $offset = 1;# another option: ($C{IR} =~ /IR/i) ? 1 : 0; foreach $pfile(@files) { print EXA "READ PRZM $pfile\n"; if($C{RUNOFF} =~ /none/i) { # insert runoff into .EXA file # EXAMS 2.98 sets monthly runoff based on the P2E-C*.d* files, so to # keep the pond with zero flow, we need to set STFLO to 0 here. print EXA "set STFLO(1,*) = 0.0\n"; print EXA "set EVAP(*,*) = 0.0\n"; print EXA "set NPSFL(*,*)=0.0\n"; print EXA "set NPSED(*,*)=0.0\n"; print EXA "set RAIN(*) = 0.0\n"; } elsif($C{RUNOFF} =~ /total/i) { # monthly runoff is default in 2.98 print EXA "SET STFLO(1,*) = $rf\n"; print EXA "set EVAP(*,*) = 0.0\n"; print EXA "set NPSFL(*,*)=0.0\n"; print EXA "set NPSED(*,*)=0.0\n"; print EXA "set RAIN(*) = 0.0\n"; } else { # assume this means monthly runoff, which EXAMS 2.98 adds print EXA "SET STFLO(1,*) = 0.0\n"; # runoff as NPSFL print EXA "set EVAP(*,*) = 0.0\n"; print EXA "set NPSED(*,*)=0.0\n"; print EXA "set RAIN(*) = 0.0\n"; } print EXA ($run++ ? "CONTINUE" : "RUN"), "\n"; $yr++; } # Temperatures used to be set here, but EXAMS 2.98 sets them, so this was # removed. (18 Sept. 2002). print EXA "QUIT"; close EXA; } ############################################### sub WriteAnex { open ANEX, ">anex.txt" or die "anex not available at this time - stay tuned\n"; print ANEX $C{IR} eq "IR" ? "1\n" : "0\n"; foreach $name("KBACW", "KBACS", "KDP") { print ANEX $C{$name} > 0 ? $C{$name} : 1.0E22, "\n"; } # For hydrolysis look in the %Hydrol hash and use the middle or the only # hydrolysis rate in it. my @k = sort keys %Hydrol; my $hyd; if(@k == 3) { $hyd = $Hydrol{$k[1]}; # use the middle hydrolysis value if there are 3 } elsif(@k == 1) { $hyd = $Hydrol{$k[0]}; # ... and the only value if there is only 1. } else { $hyd = 1.0E22; print "Error in Anex: bad number of hydrolysis values. Assuming stable\n"; } print ANEX $hyd > 0 ? $hyd : 1.0E22, "\n"; my $kps = $C{KD} ne "" ? $C{KD} * 25 : $C{KOC};# assume an OC fraction of 4% print ANEX "$kps\n"; my $rflat = 40; print ANEX "$rflat\n"; print ANEX "$C{MWT}\n"; print ANEX "$C{SOL}\n"; print ANEX "$C{VAPR}\n"; print ANEX $filepath."\\metfiles\\".$metfilename; close ANEX; print "Anex.txt written\n"; } ############################################### # Gets the total runoff for the run from the .zts file (it's just the # runoff value from the last line of the file, assuming the file prints # cumulative runoff values). Converts to m^3/hr sub GetTotalRunoff { my $flow; my $ztsfile = $oldOS ? substr($outputname,0,8) : $outputname; if(open ZTS, $ztsfile.".zts") { # PRZM truncates filename my $line; my $n = 1; # the first line will be read in the header-skip loop. while() { # skip header. header is any line with a letter after last unless(substr($_,8) =~ /[A-DF-Z\/]/i); # position 8. } while() { # go through the data portion of the file $n++; # count the lines. $line = $_; } close ZTS; my ($afield, $hl) = FieldSize(); my $convertFactor = 100.0*$afield/24.0/$n; $flow = substr($line,36,14) * $convertFactor; } else { print STDERR "Error opening $outputname.zts. Flow set to zero\n"; $flow = 0; } print "returning flow = $flow\n"; return $flow; } ############################################### # returns an array (Jan=index 0) of monthly runoff values # The P2E file gives cm/day. Before return the units are changed to m^3/hr # by dividing by 240 (unit conversion of 1/2400 and a field size as read in # from the spreadsheet file. sub GetMonthlyRunoff { my ($P2Efile) = @_; my @r = (0) x 12; my @daysInMonth = (31,28,31,30,31,30,31,31,30,31,30,31); open(PRZMFILE,$P2Efile) or die("Couldn't open $P2Efile to get runoff\n"); while() { last if(/^\s!!\*/); } # skip the top part of the file ; # skip next line $_ = ; my @temp = split(/\s+/); # The last field is the No. of apps for(my $i=0 ; $i<$temp[$#temp] ; $i++) { # skip the application data ; } while() { # read in the data my @data = split(/\s+/); shift @data unless $data[0]; # remove a blank field at the beginning, if there is one $r[$data[0]-1] += $data[2]; #add to total for month } my $conv_factor = $P{AFIELD}*100.0; for(my $i=0 ; $i<@r ; $i++) { # change units and account for field size # this converts to m^3/ hour by dividing by the total hours in the month # leap years are ignored. $r[$i] *= $conv_factor / ($daysInMonth[$i]*24); } close PRZMFILE; return @r; } ################################################################## # This subroutine reformats the FGETSEXP.XMS file sub TimeSeries { open EXMOUT, "FGETSEXP.XMS" or push @errors, "Couldn't open fgetsexp.xms"; my @one; # two arrays for holding data from the file my @two; my @results; while() { last unless(/^\s*!/); # skip the header lines (any lines starting with !) } chomp; s/^\s+//; # remove any leading whitespace in the line @one = split(/\s+/); while() { chomp; # read in a line into @two if the 1st field is the same as the 1st field # of the previous line read, copy it to @one, otherwise print out @one and # then copy @two to it. s/^\s+//; # remove any leading whitespace in the line @two = split(/\s+/); if($two[0] eq $one[0]) { @one = @two; } else { push @results, "$one[3]\t$one[4]"; # col water and benthic @one = @two; $line++; # increment the line only after adding to the column } } push @results, "$one[3]\t$one[4]"; # add in the last day of the sim. close EXMOUT; # system "del FGET*.XMS"; print "writing Time Series file $outputname_TS.out\n"; open TS, ">$outputname"."_TS.out" or push @errors, "couldn't open time series file for writing\n"; print TS "day\tdate\tconc (ppm)\tbenthic conc (ppm)\n"; my $startYr = "19".substr($startdate,-2,2); for(my $i=0 ; $i<@results ; $i++) { my $date = DateConvert($i+1,$startYr); print TS "$i\t$date\t$results[$i]\n"; } close TS; } ############################ # Two utilities used by the TimeSeries routine above for date conversion sub DateConvert { my ($daynum,$baseyear) = @_; my @daysInMonth = (31,28,31,30,31,30,31,31,30,31,30,31); my($yr,$mn,$dy) = ($baseyear,1,1); my $diy = DaysInYear($yr); while($daynum > $diy) { $daynum -= $diy; $diy = DaysInYear(++$yr); } for($mn=1 ; $mn < 12 ; $mn++) { my $dim = $daysInMonth[$mn-1]; $dim = 29 if($dim eq 28 and $diy eq 366); # leap day last if $daynum <= $dim; $daynum -= $dim; } $dy = $daynum; my $date = sprintf("%02d/%02d/%04d",$mn,$dy,$yr); return $date; } sub DaysInYear { my $y = shift; if($y % 4 eq 0) { if($y % 100 eq 0# ($,,$") = ("\t", "\t"); ) { if($y % 400 eq 0) { return 366; } else { return 365; } } return 366; } return 365; } ################################################# # Does pretty much what the Table20 program does. It extracts # data from table 20 of report.xms and writes it to report.txt sub Table20 { my $fn = $outputname.shift @_; my @sections = @_;# keywords for subsections of table 20from which to extract # data. There are main subsections in Table 20 for Direct # Exposure, Trophic Exposure and Total Media concentrations # and subsubsections in each for water columna and benthic print "writing $fn.out; benthic = $benthic\n"; open(my $REP, ">$fn.out") or die("couldn't open $fn.out for writing -$?\n"); # print some summary information about the run print $REP "stored as $fn.out\n"; print $REP "Chemical: $P{PSTNAM}\n"; my $fmodtime = fileModTime($filepath."\\PRZMenv\\".$PRZMenv); print $REP "PRZM environment: $PRZMenv\tmodified $fmodtime\n"; my $fmodtime = fileModTime($filepath."\\EXAMSenv\\".$EXAMSenv); print $REP "EXAMS environment: $EXAMSenv\tmodified $fmodtime\n"; my $fmodtime = fileModTime($filepath."\\metfiles\\".$metfilename); print $REP "Metfile: $metfilename\tmodified $fmodtime\n"; print $REP "$sections[1] segment concentrations (ppb)\n\n"; open(XMS, "report.xms") or die("couldn't open report.xms -$?\n"); my ($sep1, $sep2) = ($,,$"); # ($,,$") = ("\t", "\t"); my $i = 0; my @peak = (); my @hr96 = (); # initialize these arrays and make them local my @dy21 = (); my @dy60 = (); my @dy90 = (); my @yrly = (); my @d = (); print $REP "Year\tPeak\t96 hr\t21 Day\t60 Day\t90 Day\tYearly\n"; while() { if(/Table 20.*(\d{4})/) { $year = $1; foreach my $seclabel(@sections) { # skip to subsection while() { last if(/$seclabel/i); } } while() { chomp; chomp; if(/Mean/i) { $data = $'; last; } } while() { chomp; if(/Peak/i) { my @pk = split(/\s+/, $'); # all these values are the same substr($data,0,0) = $pk[$#pk]; # put one before the $data string @d = split(/\s+/,$data); for($i = 0 ; $i<@d ; $i++) { $d[$i] *= 1000; #convert ppm to ppb } push(@peak, $d[0]); push(@hr96, $d[1]); push(@dy21, $d[2]); push(@dy60, $d[3]); push(@dy90, $d[4]); push(@yrly, $d[5]); print $REP $year; for my $d(@d) {print $REP "\t$d";} print $REP "\n"; last; } } } } @peak = sort {$b <=> $a} @peak; @hr96 = sort {$b <=> $a} @hr96; @dy21 = sort {$b <=> $a} @dy21; @dy60 = sort {$b <=> $a} @dy60; @dy90 = sort {$b <=> $a} @dy90; @yrly = sort {$b <=> $a} @yrly; print $REP "\nSorted results\n"; print $REP "Prob.\tPeak\t96 hr\t21 Day\t60 Day\t90 Day\tYearly\n"; my $nyrs = @peak; # or any of the other arrays for($i = 0 ; $i<$nyrs ; $i++) { my $prob = ($i+1)/($nyrs+1); print $REP "$prob\t$peak[$i]\t$hr96[$i]\t$dy21[$i]\t$dy60[$i]\t$dy90[$i]\t$yrly[$i]\n"; }# ($,,$") = ("\t", "\t"); # calculate 10th %-ile value(s): my @pr = (0.1); push @pr, split(/\s/,$C{'%-ILE'}) if exists $C{'%-ILE'}; @pr = sort {$a <=> $b} @pr; print $REP "\n"; foreach(@pr) { $p = $_*($nyrs+1) - 1; # ($,,$") = ("\t", "\t"); # the "index" to interpolate to (-1 for zero offset) my($h, $l) = (int($p), int($p)+1); # indexes to interpolate between my $interpFactor = ($p-$h)/($l-$h); @interps = ($_); push @interps, $interpFactor*($peak[$l]-$peak[$h])+$peak[$h]; push @interps, $interpFactor*($hr96[$l]-$hr96[$h])+$hr96[$h]; push @interps, $interpFactor*($dy21[$l]-$dy21[$h])+$dy21[$h]; push @interps, $interpFactor*($dy60[$l]-$dy60[$h])+$dy60[$h]; push @interps, $interpFactor*($dy90[$l]-$dy90[$h])+$dy90[$h]; push @interps, $interpFactor*($yrly[$l]-$yrly[$h])+$yrly[$h]; print $REP "$interps[0]"; for(my $i=1 ; $i<@interps ; $i++) { print $REP "\t$interps[$i]";} print $REP "\n"; } my $overallAverage = 0; foreach $v(@yrly) { $overallAverage += $v; } $overallAverage /= @yrly; print $REP "\t\t\t\t\tAverage of yearly averages:\t$overallAverage\n"; # ($,,$") = ($sep1, $sep2); close XMS; print $REP "\nInputs generated by $datemessage\n"; print $REP "\nData used for this run:\n"; writeChemRunInputs($REP); close REP; }