#!/usr/bin/perl -w

#  Create Colorado groundwater input data files
#  based on well by well acreage assignment
#
#  Willem A Schreuder
#  May 22, 2006

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

#  Counties
my @cty = ('CHEYENNE','KITCARSON','LINCOLN','LOGAN','PHILLIPS','SEDGWICK','WASHINGTON','YUMA');
#  Monthly pumping/recharge distribution
my %frac = ('01'=>0,'02'=>0,'03'=>0,'04'=>0,'05'=>0,'06'=>0.147,'07'=>0.320,'08'=>0.328,'09'=>0.205,'10'=>0,'11'=>0,'12'=>0);

#
#  Model domain
#
my ($Ni,$Nj) = (165,326);
my ($X0,$Y0) = (266023,14092806);
my ($DX,$DY) = (5280,5280);
my ($X1,$Y1) = ($X0+$Nj*$DX,$Y0+$Ni*$DY);
sub off {my ($i,$j) = @_; return ($i-1)*$Nj+($j-1);}
my @ibound;
open(DAT , "<../static/02.ibound") || die "Cannot open file ../static/02.ibound\n";
foreach my $i (1 .. $Ni)
{
   my $line = <DAT>; chomp $line;
   foreach my $j (1 .. $Nj)
   {
      $ibound[off($i,$j)] = substr($line,2*($j-1),2);
   }
}
close(DAT);

#
#  Get line of tagged data and scale
#
sub get
{
   my ($key,$n,$div) = @_;
   my $line = <DAT>;
   chomp $line;
   my ($tag,@line) = split(' ',$line);
   ($n==@line) || die "Incorrect item count reading $key\n";
   ($key eq uc($tag)) || die "Expected $key, found $tag\n";
   return map {$_/$div} @line;
}

#  Usage
my %args = GetOpt('f:d:',f=>undef,d=>undef);
(@ARGV>0) || die "Usage: $0 <year>\n";
foreach my $year (@ARGV)
{
   #
   #  Summary
   #
   my %sum;
   foreach my $cty (@cty,'SPECIAL','CCP')
   {
      $sum{$cty} = {PUMP=>0,RECH=>0,AREA=>0};
   }
   #
   #  Get CIR and efficiency data
   #
   my (%CIR,%pump,%rech);
   my $file = $args{f} ? $args{f} : "../data/co/$year.gw";
   my $dir  = $args{d} ? $args{d} : "co";
   open(DAT , "<$file") || die "Cannot open file $file\n";
   my ($deficit)             = get('DEFICIT'  ,1,100);
   my ($FarmEFFspr,$PERCspr) = get('SPRINKLER',2,100);
   my ($FarmEFFfld,$PERCfld) = get('FLOOD'    ,2,100);
   @CIR{@cty} = get('NETCIR',scalar(@cty),12); 
   #  Pumping by county
   foreach my $cty (@cty)
   {
      $pump{SPR}{$cty} = $deficit * $CIR{$cty} / $FarmEFFspr;
      $pump{FLD}{$cty} = $deficit * $CIR{$cty} / $FarmEFFfld;
      $rech{SPR}{$cty} = $PERCspr*$pump{SPR}{$cty};
      $rech{FLD}{$cty} = $PERCfld*$pump{FLD}{$cty};
   }

   #
   #  Read well locations and acreage
   #  Calculate pumping, recharge and area from well by well values
   #
   my @agw = (0) x ($Ni*$Nj);
   my (%pmp,%rcg);
   foreach my $mo (sort keys %frac)
   {
      @{$pmp{$mo}} = (0) x ($Ni*$Nj);
      @{$rcg{$mo}} = (0) x ($Ni*$Nj);
   }
   while (my $line = <DAT>)
   {
      #  Skip comment lines
      ($line =~ /^#/) && next;
      #  Read well location and acres
      chomp $line;
      my ($x,$y,$cty,$Afld,$Aspr,$Area) = split(' ',$line);
      #  Discard wells outside model domain
      ($x<$X0 || $x>$X1 || $y<$Y0 || $y>$Y1) && next;
      #  Discard wells in dead cells
      my ($i,$j) = ($Ni-int(($y-$Y0)/$DY),int(($x-$X0)/$DX)+1);
      my $off = off($i,$j);
      ($ibound[$off] == 0 && $cty eq '*') && die "Transfer from dead cell\n";
      ($ibound[$off] == 0) && next;
      #  CCP has monthly pumping
      if ($cty eq 'CCP')
      {
         my (undef,undef,undef,@pump) = split(' ',$line);
         foreach my $mo (sort keys %frac)
         {
            my $pump = shift @pump;
            $sum{CCP}{PUMP} += $pump;
            $pmp{$mo}[$off] += $pump;;
         }
      }
      # Others have annual pumping
      else
      {
         my ($area,$pump,$rech);
         #  Special case
         if ($cty eq '*')
         {
            $cty = 'SPECIAL';
            $pump = $Afld;  #  Flood field is the total pumping (af)
            $rech = $Aspr;  #  Sprinkler field is the returns (af)
            $area = $Area;  #  Appropriated field as actual acres
         }
         #  Irrigation well by county
         else
         {
            #  Check for valid county
            $cty = uc($cty);
            exists($CIR{$cty}) || die "Unknown county $cty\n";
            #  Pumping and return flow
            $area = $Aspr + $Afld;
            $pump = $pump{SPR}{$cty}*$Aspr + $pump{FLD}{$cty}*$Afld;
            $rech = $rech{SPR}{$cty}*$Aspr + $rech{FLD}{$cty}*$Afld;
         }
         #  Assign to cells by month
         $agw[$off] += $area;
         while (my ($mo,$frac) = each %frac)
         {
            $pmp{$mo}[$off] += $frac*$pump;
            $rcg{$mo}[$off] += $frac*$rech;
         }
         #  Summarize values by county
         $sum{$cty}{PUMP} += $pump;
         $sum{$cty}{RECH} += $rech;
         $sum{$cty}{AREA} += $Aspr + $Afld;
      }
   }
   close(DAT);

   #
   #  Save to file
   #
   sub save
   {
      my ($file,@x) = @_;
      #  Skip empty files
      my $sum=0;
      foreach my $x (@x)
      {
         $sum += abs($x);
      }
      $sum || return;
      #  Dump values to file
      open(DAT , ">$file") || die "Cannot open file $file\n";
      print DAT join("\n",@x)."\n";
      close(DAT);
   }

   #  Annual acres
   save("$dir/$year.agw",@agw);
   #  Monthly pumping and recharge
   while (my ($mo,$frac) = each %frac)
   {
      save("$dir/$year.$mo.pmp",@{$pmp{$mo}});
      save("$dir/$year.$mo.rcg",@{$rcg{$mo}});
   }

   #  Total
   foreach my $cty (@cty,'SPECIAL','CCP')
   {
      foreach my $var ('PUMP','RECH','AREA')
      {
         $sum{TOTAL}{$var} += $sum{$cty}{$var};
      }
   }
   #  Print summary
   printf "%-12s %8s %8s %8s\n" , 'County' , 'Pumping' , 'Returns' , 'Acres';
   foreach my $cty (@cty,'SPECIAL','CCP','TOTAL')
   {
      printf "%-12s %8d %8d %8d\n" , $cty , @{$sum{$cty}}{'PUMP','RECH','AREA'};
   }
}
