#!/usr/bin/perl -w

#  This program takes the Nebraska groundwater and surface water data,
#  both exclusive and comingled, and builds data sets for the preprocessor.
#  Well data is converted into a format suitable for mkgw to generate
#  pumping and groundwater returns.
#
#  Private river canals were previously handled by two files,
#  PrivateCanalLosses_to_Fields_NE_061003.dbf and SWexLand_Revised.dbf. These
#  canals are now handled in the surface water file.  In the surface water file,
#  the column "Volume" for private canals represents the diversion. Return flow
#  from the private canals is specified system and is always 0.4 times the
#  diversion as specified in the final model report and agreed by the RRCA
#  committee.  "Volume" for all other canals represents the volume of
#  application to the fields.
#
#  In addition the River Pumpers file, RiverPumpers_NE_061003.dbf, was
#  eliminated from use.  This is now handled by the surface water comingled and
#  exclusive files using a system number greater than 9000.
#
#  This file was revised by Chuck Spalding from the original MKNEDAT as
#  contained on the July 1,2003 DVD-- Appendix A of the Republican River
#  Compact Admistration Groundwater model June 30 2003.  The original MKNEDAT
#  file was created on 6/29/03 @10:01 PM.
#
#  Further revisions to generalize program for all years from 2000 by Willem
#  Schreuder 12/30/2003
#
#  Modified 4/24/2014 to include augmentation pumping (Willem Schreuder)

use strict;
use lib '/pm/grflib/perl';
use GetOpt;
use DBF;

my %args = GetOpt('d:',d=>'../data/ne/export');
my $dir = $args{d};
(@ARGV>0) || die "Usage: $0 [-d <dir>] <year(s)>\n";

#  MODFLOW domain
my $Nrow = 165;
my $Ncol = 326;
my $N    = $Nrow*$Ncol;

#
#  Irrigation distribution
#  Maps annual volume to months
#
my %irrdist = (6=>0.147,7=>0.320,8=>0.328,9=>0.205);
my %month = (1=>'Jan',2=>'Feb',3=>'Mar',4=>'Apr',5=>'May',6=>'Jun',7=>'Jul',8=>'Aug',9=>'Sep',10=>'Oct',11=>'Nov',12=>'Dec');

#
#  River Pumpers return flow fractions
#
my $RPrf = 0.25;

#
#  Groundwater pumping return flow fraction 0.2 from 2000
#
my $GWrf = 0.2;

#
#  SW system return flow fractions
#
my %SWrf = (
   #-----------Federal----------------
   800=>0.246,   #  Culbertson
   801=>0.246,   #  Culbertson Extension
   802=>0.246,   #  Meeker-Driftwood
   803=>0.246,   #  Bartley
   804=>0.246,   #  Red Willow
   805=>0.246,   #  Cambridge
   807=>0.287,   #  Franklin
   808=>0.287,   #  Naponee
   809=>0.287,   #  Franklin Pump
   810=>0.2542,  #  Superior
   811=>0.1886,  #  Courtland
   #-----------Private----------------
   812=>0.4,     #  Haigler
   813=>0.4,     #  Champion
   814=>0.4,     #  Riverside
   815=>0.4,     #  Unknown - new in 2008
   #-----------CNPPID----------------
   816=>0.4,     #  Jeffrey Hydro
   817=>0.4,     #  Jeffrey Hydro to J1 Hydro
   818=>0.4,     #  J2 to the Platte
   819=>0.3,     #  Phelps |
   820=>0.3,     #  E65    | Efficiency after 2000
   821=>0.3,     #  E67    |
   #-----Platte Private Wells--------
   822=>0.4,     #  Sixmile
   823=>0.4,     #  Thirtymile
   824=>0.4,     #  Orchard
   826=>0.4,     #  Western
   #-------------NPPD----------------
   825=>0.4,     #  NPPD
   #-----Miscellaneous-----
   899=>0.4      #  Private Canal Losses to Fields
   );

#
#  Offset from row and column
#  Store FORTRAN style (column major order)
#
sub offset
{
   my ($row,$col) = @_;
   return $Ncol*($row-1)+($col-1);
}

#
#  Store volumes for rrpp unless all values are zero
#
sub store
{
   my ($yr,$mo,$ext,@vol) = @_;
   defined($mo) && ($mo = sprintf '%.2d' , $mo);
   #  Count nonzero values
   my $n=0;
   foreach my $vol (@vol)
   {
     ($vol != 0) && $n++;
   }
   ($n>0) || return;
   #  Save file
   my $file = defined($mo) ? "ne/$yr.$mo.$ext" : "ne/$yr.$ext";
   open(DAT , ">$file") || die "Cannot open file $file\n";
   foreach my $vol (@vol)
   {
      print DAT "$vol\n";
   }
   close(DAT);
}

#
#  Open all the databases and index on year
#
my %db;
$db{GW}{EX} = DBF->open("$dir/GWOnlyExport.dbf",'DIRECT'); $db{GW}{EX}->index('Year','N');
$db{GW}{CO} = DBF->open("$dir/GWCoExport.dbf"  ,'DIRECT'); $db{GW}{CO}->index('Year','N');
$db{SW}{EX} = DBF->open("$dir/SWOnlyExport.dbf",'DIRECT'); $db{SW}{EX}->index('Year','N');
$db{SW}{CO} = DBF->open("$dir/SWCoExport.dbf"  ,'DIRECT'); $db{SW}{CO}->index('Year','N');
$db{AW}     = DBF->open("$dir/AugExport.dbf"   ,'DIRECT'); $db{AW}->index('Year','N');

#
#  Process data for requested years
#
foreach my $yr (@ARGV)
{
   my (%cty,%volm,%area,%rate);
   #  Process all the wells
   foreach my $typ ('EX','CO')
   {
      #  Get well data for this year - need county totals only
      foreach my $rec ($db{GW}{$typ}->select('Year',$yr,'N','ROW','COLUMN','County','Area','Volume'))
      {
         my ($row,$col,$cty,$acr,$pmp) = @{$rec}{'ROW','COLUMN','County','Area','Volume'};
         #  Remember county
         $cty{$cty} = 1;
         #  Accumulate pumping and acres by type, county and year
         $volm{GW}{$typ}{$cty} += $pmp;
         $area{GW}{$typ}{$cty} += $acr;
      }
   }
   #  Process surface acres - need county totals only
   foreach my $typ ('EX','CO')
   {
      #  Get surface data for this year
      foreach my $rec ($db{SW}{$typ}->select('Year',$yr,'N','County','Area','Volume'))
      {
         my ($cty,$acr,$div) = @{$rec}{'County','Area','Volume'};
         #  Remember county
         $cty{$cty} = 1;
         #  Accumulate diversions and acres by type, county and year
         $volm{SW}{$typ}{$cty} += $div;
         $area{SW}{$typ}{$cty} += $acr;
      }
   }

   #  Make sorted county list
   my @cty = sort keys %cty;

   #
   #  Set missing values to zero
   #  Check to make sure GW and SW data agree on comingled acres
   #  Translate pumping into rates for comingled
   #
   foreach my $cty (@cty)
   {
      #  Set missing years to zero
      foreach my $src ('GW','SW')
      {
         foreach my $typ ('EX','CO')
         {
            exists($volm{$src}{$typ}{$cty}) || ($volm{$src}{$typ}{$cty} = 0);
            exists($area{$src}{$typ}{$cty}) || ($area{$src}{$typ}{$cty} = 0);
         }
      }
      #  Check that SW & GW comingled areas match to withing 10 acre-feet
      (abs($area{SW}{CO}{$cty}-$area{GW}{CO}{$cty})>10) && die "Error SW-GW comingling area mismatch $cty $area{SW}{CO}{$cty} $area{GW}{CO}{$cty}\n";
      #  Calculate rate of application - 0 if none
      $rate{GW}{CO}{$cty} = ($area{GW}{CO}{$cty}>0) ? $volm{GW}{CO}{$cty}/$area{GW}{CO}{$cty} : 0;
      $rate{GW}{EX}{$cty} = ($area{GW}{EX}{$cty}>0) ? $volm{GW}{EX}{$cty}/$area{GW}{EX}{$cty} : 0;
      $rate{SW}{CO}{$cty} = ($area{SW}{CO}{$cty}>0) ? $volm{SW}{CO}{$cty}/$area{SW}{CO}{$cty} : 0;
      $rate{SW}{EX}{$cty} = ($area{SW}{EX}{$cty}>0) ? $volm{SW}{EX}{$cty}/$area{SW}{EX}{$cty} : 0;
   }

   #
   #  Save county summary
   #
   open(DAT , ">necounty.dat") || die "Cannot open file necounty.dat\n";
   printf DAT '%-12s' , $yr;
   printf DAT '%12s' , 'Vol GW Ex';
   printf DAT '%12s' , 'Vol GW Co';
   printf DAT '%12s' , 'Vol SW Ex';
   printf DAT '%12s' , 'Vol SW Co';
   printf DAT '%12s' , 'Rate GW Ex';
   printf DAT '%12s' , 'Rate GW Co';
   printf DAT '%12s' , 'Rate SW Ex';
   printf DAT '%12s' , 'Rate SW Co';
   printf DAT '%12s' , 'Area GW Ex';
   printf DAT '%12s' , 'Area SW Ex';
   printf DAT '%12s' , 'Area SW Co';
   print DAT "\n";
   foreach my $cty (@cty)
   {
      printf DAT "%-12s" , $cty;
      printf DAT "%12.2f" ,$volm{GW}{EX}{$cty};
      printf DAT "%12.2f" ,$volm{GW}{CO}{$cty};
      printf DAT "%12.2f" ,$volm{SW}{EX}{$cty};
      printf DAT "%12.2f" ,$volm{SW}{CO}{$cty};
      printf DAT "%12.2f" ,$rate{GW}{EX}{$cty};
      printf DAT "%12.2f" ,$rate{GW}{CO}{$cty};
      printf DAT "%12.2f" ,$rate{SW}{EX}{$cty};
      printf DAT "%12.2f" ,$rate{SW}{CO}{$cty};
      printf DAT "%12.2f" ,$area{GW}{EX}{$cty};
      printf DAT "%12.2f" ,$area{SW}{EX}{$cty};
      printf DAT "%12.2f" ,$area{SW}{CO}{$cty};
      print DAT "\n";
   }

   #
   #  Calculate cell by cell pumping, recharge, area and rate
   #
   #  Set volumes amd areas to 0
   my (%acr,%ret,%pmp);
   @{$pmp{Ann}} = (0) x ($N);
   @{$acr{GW}} = @{$acr{SW}} = @{$acr{CO}} = (0) x ($N);
   @{$ret{GW}} = @{$ret{SW}} = (0) x ($N);
   #  Get all the surface water diversions
   foreach my $typ ('EX','CO')
   {
      #  Process individual lands for this year
      foreach my $rec ($db{SW}{$typ}->select('Year',$yr,'N','ROW','COLUMN','County','Area','Volume','System'))
      {
         my ($row,$col,$cty,$acr,$div,$sys) = @{$rec}{'ROW','COLUMN','County','Area','Volume','System'};
         my $l = offset($row,$col);
         #  For comingled lands, add groundwater to surface water
         #  GW = comingled rate for county and year times area
         my $pmp = ($typ eq 'CO') ? $rate{GW}{CO}{$cty}*$acr : 0;
         #  No water this year - skip the rest
         ($div+$pmp>0) || next;
         #  River pumpers have system codes > 999
         ($sys>999 || exists($SWrf{$sys})) || die "Unknown SW system $sys\n";
         my $SWrf = ($sys>999) ? $RPrf : $SWrf{$sys};
         #  Add area to SW or CO lands
         $acr{($typ eq 'EX') ? 'SW' : 'CO'}[$l] += $acr;
         #  Add return flows to GW and SW volumes
         $ret{GW}[$l] += $SWrf*$pmp;
         $ret{SW}[$l] += $SWrf*$div;
      }
   }
   #  Get all the wells - static monthly distribution
   foreach my $typ ('EX','CO')
   {
      #  Get well data for this year
      foreach my $rec ($db{GW}{$typ}->select('Year',$yr,'N','ROW','COLUMN','County','Area','Volume'))
      {
         my ($row,$col,$cty,$acr,$pmp) = @{$rec}{'ROW','COLUMN','County','Area','Volume'};
         my $l = offset($row,$col);
         #  No water this year - skip the rest
         #  After 2014 zero pumping is culled and this represents rollover acreage
         ($yr>=2014 || $pmp>0) || next;
         #  Add pumping to cell
         $pmp{Ann}[$l] += $pmp;
         #  Add area, volume and return flows to groundwater exclusive lands
         if ($typ eq 'EX')
         {
            $acr{GW}[$l] += $acr;       #  Record GW acres
            $ret{GW}[$l] += $pmp*$GWrf; #  Record GW returns
         }
      }
   }
   #  Augmentation pumping by month
   my %aug;
   foreach my $mon (values %month)
   {
      @{$pmp{$mon}} = (0) x ($N);
      foreach my $rec ($db{AW}->select('Year',$yr,'N','ROW','COLUMN','Project',$mon))
      {
         my ($row,$col,$prj,$pmp) = @{$rec}{'ROW','COLUMN','Project',$mon};
         my $l = offset($row,$col);
         #  Add pumping to cell
         $pmp{$mon}[$l] += $pmp;
         #  Summarize augmentation
         $aug{$prj}{$mon} += $pmp;
      }
   }
   #  Save augmenttion pumping summary
   if (%aug)
   {
      open(DAT , '>augmentation.dat') || die "Cannot open file augmentation.dat\n";
      printf DAT "%-12s" , $yr;
      foreach my $mo (1 .. 12)
      {
         printf DAT "%10s" , $month{$mo};
      }
      printf DAT "%10s\n" , 'Total';
      foreach my $prj (sort keys %aug)
      {
         printf DAT "%-12s" , $prj;
         my $sum = 0;
         foreach my $mo (1 .. 12)
         {
            $sum += $aug{$prj}{$month{$mo}};
            printf DAT "%10.2f" , $aug{$prj}{$month{$mo}};
         }
         printf DAT "%10.2f\n" , $sum;
      }
      close(DAT);
   }
   #  Save annual irrigated area files
   foreach my $typ ('SW','GW','CO')
   {
      store($yr,undef,'a'.lc($typ),@{$acr{$typ}});
   }
   #  Save monthly pumping and returns
   foreach my $mo (1 .. 12)
   {
      my $mon = $month{$mo};
      #  Agricultural pumping and returns
      if (exists($irrdist{$mo}))
      {
         #  Add ag pumping to augmentation pumping
         for (my $l=0;$l<$N;$l++)
         {
            $pmp{$mon}[$l] += $irrdist{$mo}*$pmp{Ann}[$l];
         }
         #  Save returns
         store($yr,$mo,'rcg',map {$irrdist{$mo}*$_} @{$ret{GW}});
         store($yr,$mo,'rcs',map {$irrdist{$mo}*$_} @{$ret{SW}});
      }
      #  Total pumping
      store($yr,$mo,'pmp',@{$pmp{$mon}});
   }
}
